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LIST OF SYMBOLS 


In the notation associated with this report, a number of 
principles have been adhered to. 

(a) The symbols A'., B, C, H, E are usually matrix quantities 
of dimensions (n x m) or (n x n). 

(b) The symbols i, j, k, 1, m, n are usually constants. 

(c) The symbols U, V., W, X, Y, Z. are usually vector 
quantities of dimensions (n x l). 

(d) Lower case symbols are single quantities, while upper case 
symbols usually represent arrays. 

(e) Where possible, the normally-accepted nomenclature is used 
to represent the quantity. 

(f) Where a symbol has more than one meaning, the meaning o 

assigned to a particular section is indicated. 

************ 

the "calibrated focal length of the camera 
x and y coordinates of an image point 

x and y coordinates of the principal point 

image displacements 




x and y image velocities 


x and y instantaneous image velocities 
the exposure epoch 

the attitude of the camera axis. The order of the 
rotations is tertiary, secondary, primary, respectively 


the rates of change of camera axis attitude 


a vector of coordinates representing the photostation 

lx* 


X o = 


a vector of object space coordinates 

fx' 


X = 

p 


a vector of velocity components of the photostation 

fv 

(used in Chapter 2) 


V = 


V* 


the general three-dimensional rotation matrix 


< . 


M = R^R^E^ (used in Chapter 2, Appendices A, F) • 



M. 

x 


the ith row of the general rotation matrix M: 

i » 1, 2, 3 


the augmented rotation matrix Used in Chapter 2, 


N = R*Rv. . , RiRx . R, R.- , 
* fra t 9 9*t w <*>nt 


Appendices A, F 


the ith row of the augmented rotation matrix N: 
i - If 2, 3 


the derivative of the ith row of N with respect to <=kt 


a vector norm 


a matrix norm 


the Turing M number 


(used in Chapter 3) 


the Turing N number 


(used in Chapter 3) 


the Hadamard Condition number 


(used in Chapter 3) 


the Todd P number 


(used in Chapter 3) 


the weight matrix P = £ 


(used in Chapter 4) 


the approximate values of the parameters 


the adjusted values of the parameters 


the adjusted values of the observations 


the observed values of the observations 



WW wffl <] W4 w<j 


approximate values of the observations computed 
through the model 


a vector of residuals V = L. - L (used in Chapter 4) 

b a 

a discrepancy vector £ = - L q 


a vector of Lagrange multipliers 


partial differential matrix of function with respect 
to observed quantities. 


partial differential matrix of function with respect 
to parameters 


partial differential matrix of function with respect 
to elements of exterior orientation 


partial differential matrix of function with respect 
to survey parameters 


the correction vector to be applied to the approximate 

values to obtain the adjusted values X = X + A 

u a oo 


the correction vector for the elements of exterior 


orientation 


the correction vector for the survey parameters 


the variance-covariance matrix 


x 



the ij element of the variance-covariance matrix 


the general symbolic form of the reduced matrix 
system 



1 . INTRODUCTION 


1.0 Scope and Objectives 

Two principal types of camera are currently available to photo- 
grammetrists . The first type is, from a historical viewpoint, the 
ideal camera for a mapping photogrammetrist,. since it usually com- 
bines low objective distortion characteristics with a high degree 
of stability. That is, its laboratory calibration can be assumed 
to apply under a wide range of conditions for an extended time 
period. The classical cameras of this type are the Wild RC and 
the Zeiss RMK series of aerial cameras. 

The second type of camera is comparatively new to photogram- 
metry. It is characterized by high image quality and information 
content, but not so stable geometry'. - .The class is typified by 
the panoramic and focal plane shutter" cameras manufactured by 
Itek and Hycon. 

Unfortunately, current technology in the design and manufacture 
of large-format camera systems is such that the two types of camera 
are almost mutually exclusive. Consequently, while the first group 
of cameras has been used almost exclusively for mapping purposes, 
the second group has been used almost exclusively for photointer- 
pretation tasks. Recently, this distinction has begun to break 
down dramatically, with cameras from the second or non-metric group 
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being -used more and more in small-scale metric tasks where the 
amount of information content is paramount. 

One of the first breaks in class distinction and use occurred 
in the mid-sixties, when NASA chose a non-metric photo subsystem 
to accomplish the metric task of mapping the Apollo zone of the 
moon. Unfortunately, at that time the geometric fidelity and 
stability of non-metric cameras were incompletely understood. This 
resulted in data mis-matches and discrepancies for the lunar tri- 
angulation. 

The general trend of previous work on non-metric cameras had 
been directed towards a complete understanding of the optical qual- 
ity of the resulting image and those factors' which affected it. 
There was now a need to also understand the geometric fidelity of 
the image. This study therefore seeks to provide meaningful 
answers to the following general problems which arise when non- 
metric cameras are used for metric work. 

(a) What is the effect of image motion and image motion compen- 
sation on the location of the principal points 

(b) Is it possible to determine the corrections to the cali- 
brated values of the coordinates (defining the location of the 
principal point or the fiducial centre) in order to correct for the 
image motion compensation? 

(c) What is the effect of the focal plane shutter on the dis- 
tortion and interior geometry (interior orientation parameters)? 

(d) Can a lack of calibration information be overcome by dynam- 
ic calibration procedures incorporated in the photographic mission? 
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In view of the relevance of these studies to the particular 
needs of NASA, the Lunar Orbiter Photographic Subsystem will be 
used to typify a non -metric system. 

Chapter 2 details a new mathematical model which takes into 
consideration image motion and movement of the camera during 
the exposure interval. The model offers a rigorous relationship 
between the object space and the image space by incorporating into 
the well-known collinearity conditions velocity' terms ' for the ex- 
posure station. The model is considered general since it is pos- 
sible to apply simplifying assumptions to obtain the presently- 
accepted explanation for image motion. Unfortunately, these gains 
in generality have been accompanied by increased complexity of the 
model and consequently possible computational instability problems. 

Chapter 3 is devoted to some of the problems faced in the solu- 
tion of large systems of equations. The two principal problems 
treated are computational stability and the economic solution of 
such equations by non-iterative methods. The segment on computa- 
tional stability introduces the concept of matrix and vector norms. 
The concept is then developed. through perturbation theory into an 
acceptable method for gauging the stability of a particular solu- 
tion. The final half of the chapter investigates a number of pos- 
sible alternative methods for the solution of large matrix systems. 

Chapter k is concerned with numerical testing of the model pro- 
posed in Chapter 2. The tests are all made with simulated lunar 
data for which true values of all parameters are known. The chap- 
ter has three main subdivisions. The first details the method of 



generating the simulated data. The second section discusses the 
single photo resection tests, while the third section discusses the 
general intersection problem. Particular attention is given in all 
tests to the problem of computational stability and to the control 
of this problem by the choice of the correct weight matrix of the 
observed quantities. 

Finally Chapter 5 reviews and summarizes the work of Chapters • 
2, 3 and k in terms of the four above-mentioned problems. 

1.1 Historical Review 

The Lunar Orbiter program was conceived primarily for the ac- 
quisition of high-quality photographic data from those portions of 
the lunar surface under consideration for the Apollo program. 

Prior to the launching of Lunar Orbiter I, man's photographic 
coverage .of the moon consisted of earth-based photography by tele- 
scopes whose maximum resolution on the lunar surface was approxi- 
mately 20 metres' 1 *. A detailed discussion on the limiting factors 
for this type of coverage is given in the ACIC report ^'Department 
of Defense 1966 Selenodetic Control Data"'**. A. very small amount 
of additional coverage at greater resolution, but with very high 
geometric distortion, was obtained from the- Ranger missions. The 

reduction of this material was accomplished almost exclusively by 

2 

photometric techniques » 

Table 1 summarizes the Lunar Orbiter missions^. 



Table 1 


Summary of Lunar Orbiter Missions 


Lunar 

Orbiter 

Orbit 

Period 

(hr) 

Inclination to 
Equatorial Axis 
(deg) 

Perilune 

Altitude 

(km) 

Apolune 

Altitude 

(km) 

Date of 
Photography 

Quality of 
Photography 
Medium High 

Resolution .Resolution 

I 

3-5 

12 and 21 

60 

1850 

August 

Acceptable 

Unusable 






18-29, 1966 


IMC failed 

II 

3-5 

12 and 21 

60 

1850 

November 

Good 

Good 






18-25, 1966 



III 

3.5 

12 and 21 

60 

1850 

February 

Excellent 

Excellent 






15-25, 1967 



IV 

12 

85 

2700 

6100 

May 

Mediocre 

Excellent 






11 - 26 , 1967 

to poor 


v* 

8 . if 

85 

200 

6000 

August 

Excellent 

Excellent 


8.3 

85 

. 100 

6000 

8 - 18 , 1967 

Excellent 

Excellent 


3.2 

85 

100 

1500 


Excellent 

Excellent 


♦ 


The three values represent the three conditions under which photography was obtained 
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The success of Lunar Orbiters II and III in acquiring the prime 
Apollo landing site data, coupled with the successful Surveyor 
flights to the moon, caused a reevaluation of the Lunar Orbit er IV 
missibn. Thus the task of extrapolating the point information of 
the Surveyor missions and the forthcoming Apollo missions was to be 
done with the aid of a full lunar coverage obtained by Lunar Orbit- 
er IV. 

1.2 The Lunar Orbiter Photographic Subsystem 

The Lunar Orbiter Photographic Subsystem has been completely 

4 

described by the Eastman Kodak Company , subcontractors to The 
Boeing Company for the Lunar Orbiter program. In addition, specif- 
ic parts of the Lunar Orbiter Photographic Subsystem have heen de- 
scribed and discussed in the scientific literature, e.g. Konecny'*, 

£ rp ' 

Kosofsky , and Norman. Despite this wealth of information, it is 
appropriate to discuss here certain aspects of the system so that 
a clear understanding of the research discussed in this report is 
obtained. 

The Lunar Orbiter Photographic Subsystem incorporated two op- 
tical systems which imaged on the same 70 mm film. The 80 mm focal 
length system, often referred to as the medium resolution system, 
used a between-the-lens shutter system and a 55 x 63 mm format. 
Image Motion Compensation (INC) was possible in a single direction, 
corresponding to the direction of movement of the focal plane shut- 
ter in the 610 mm system. The magnitude of this compensation was 
determined by a V/H sensor using an unimaged segment of the high 
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resolution imagery. A mechanical reduction linkage between the 
platen of the 6l0 mm system, where the determined IMG/ was.' applied, 

T S 

and the platen of the 80 mm system provided IMG. f.or the 80 mm 
system. 

The interior orientation for the 80 mm system was determined 
by laboratory calibration techniques. The subcontractor for this 

g 

work was the Fairchild Camera Company . The calibration included 
radial and tangential distortion components at discrete points 
along each of the principal diagonals. 

The 610 mm system, often referred to as the high resolution 
system, used a focal plane shutter which traversed the short dimen- 
sion of the 55 x 219 mm format. As previously indicated, IMC 
could be applied only in the direction of motion of shutter. This 
resulted in the IMC device remaining unactivated for the Lunar 
Orbiter IV mission. Figure 1 illustrates the different configura- 
tions controlling the applicability of IMC. 

5 

The camera was only partially calibrated by Brown due to , 
operational restrictions • Thus all missions were flown with non- 
recoverable photographic systems in which both the decentering 
distortion of the photographic system and one of the coordinates 
of the principal point (y) were unknown. 

All photography derived from the Lunar Arbiter spacecraft was 
processed on board the spacecraft by the Kodak Bimat process. The 
resulting image was electronically scanned and transmitted back to 
earth where it was reassembled. Unfortunately, the framelets suf- 
fered distortion which can only be minimized. The reduction and 




SHUTTER DIRECTION 


FIGURE 1 


DIAGRAM ILLUSTRATING THE OR8IT AND CAMERA 
ORIENTATION TOR MISSIONS 1,11, III, AND IV 
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elimination of this type of distortion is a comparatively new 
problem in photo grammetry . Fortunately, a considerable effort is 
being made into understanding the causes, effects and elimination 
of this distortion. Interested readers are referred to the work 
of Wong 10 . 

Throughout this report, it is assumed that the error sources 
of this problem have been identified and therefore need not be 
considered. 



2. IMAGE MOTION AND IMAGE MOTION COMPENSATION 


2.0 Introduction 

In the late fifties and early sixties, considerable technical 
advances were made in all classes and types of cameras. In par- 
ticular a new type of camera, the reconnaissance camera, became 
firmly established. Since this camera stresses image quality over 
geometric fidelity, it is the natural complement of the survey 
camera which is geometrically stable but has relatively poor reso- 
lution. The initial need for high image-quality systems was gen- 
erated by military surveillance. Today, however, civilian needs 
for such a system are also very great, especially if the possi- 
bilities of extraterrestrial photogrammetry and photointerpreta- 
tion are to be maximized. 

Unfortunately, the operational requirements which necessitate 
these high-resolution systems are not favorable to geometric 
photogrammetry, since the photostatidn cannot be - considered as a 
unique point in space. This movement of .the photostation causes 
the image to blur. The amount of blur represents the magnitude 
of image motion present during the exposure interval. This degra- 
dation of the image cannot be tolerated and hence Image Motion 
Compensation (IMC) has been incorporated into these systems. The 
first detailed study of the degrading effects of image motion on 


10 
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a photographic image was reported in 1955 by Wolfe and Lamberts 
There are three main methods of accomplishing .IMC . They are 
as follows: 

(a) Movement of the platen-film assembly. 

(b) Movement of the lens cone. 

(c) Use of the focal plane shutter. 

Often method (c) is combined with (a) or (b). Application of 
method (a) or (b) destroys the interior geometry of the camera 
as currently conceived, while method (c) represents an infinite 
number of photographs joined side by side from a continuously 
changing photostation. For these reasons, systems which employ 
IMC have not been favored by mapping photogrammetrists. An ex- 
ception is the recently developed Fairchild KC-6A camera. 

The initial researchers in this new field of image degradation 

11 12 13 

due to motion (Wolfe and Lamberts , Trott , Eosenau ^ ) v. r ere all 

primarily concerned with the quality of the image without regard 

l4 

to its geometric position. In 19^5 > Kawachi an’d Weinflash 

derived expressions for image velocities as functions of image 

position, focal length and rotational velocities. The expressions 

were not general but provided an acceptable method of determining 

expected blur distances and hence degradation under certain given 

conditions. The work of Kawachi and Weinflash was subsequently 

15 

published in Photogrammetric Engineering . 

This article was preceded by a discussion of image motion 
resulting from translations only"^. Subsequent to these two 
papers by Kawachi, there are no published investigations on the 
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problems of image motion. Unfortunately, this could be due to 

the following statement by Kawachi^ « 

... An analytical derivation involving matrices was also 
considered since a rotation matrix transforms the coordi- 
nates of a point in one system to its coordinates in a . 
rotated coordinate system. It appears logical, there- 
fore, that matrices can be applied to obtain the image 
velocity, since the aircraft motions under consideration 
are rotational. However the matrix approach is not appro- 
priate for two reasons: (1) matrices do not describe the 
actual movement of the image point (unless additional 
translational terms are introduced to account for the 
changes in the radius vector), because the image remains 
in the film plane and hence its motion is not equivalent 
to holding the image point fixed and rotating the coor- 
dinate system; and (2) the matrix approach is not simpler 
than the geometrical approach but actually involves more 
equations. Multiplications of the matrices would show 
this . . . 

A similar statement is to be found in the Fairchild Technical Memo 
Note referenced above. 

Unfortunately, the engineers and scientists in charge of the 
Lunar Orbiter Missions chose a photographic system of the recon- 
naissance type for their mapping program. This decision was, in 
part, influenced by the resolution of such systems. However, it 
was to lead to a number of awkward problems. In particular, the 
problem of image motion and its compensation has severely limited 
the accuracy of the subsequent triangulations * since Kawachi's for- 
mulae do not interconnect the object and image spaces. These for- 
mulae are therefore not applicable for direct incorporation into 
the mathematical relationships between the two spaces. 

After failing to adequately incorporate Kawachi's formulae 
into the mathematical model, it became apparent that a new approach 
was needed. It was thus decided to investigate the feasibility of 



13 


using the collinearity conditions, since these conditions describe 
the existing relationships between object and image space systems. 


2.1 The Collinearity Condition Approach 

The well-known collinearity condition is, briefly speaking, a 

functional relationship between points in the object space and the 

image space systems. A complete derivation of this relationship 

may be found in most texts on photo gramme try, e.g. Manual of 

17 

Photogrammetry . The relationship is often written in a form 
that is pertinent to the user's particular purpose. In this 
work, it was desirable to express the collinearity condition in 
the following form: 

x - x = -f Fm. (X - X ) / M,(X - X )1 

p o tip o ' 3 p o J 


( 2 . 1 ) 


y p “ 






M_ (X - 
2 P 


X ) / M (X - X 
o ' 3 p o 


>1 


That is, given the position and attitude of the camera, it is pos- 
sible to compute for every object space point a corresponding 
image space point. There is no restriction as to what type of sur- 
face the object space field must satisfy. 


2.11 Augmentation of the Collinearity Condition Relationship 

Consider a moving photographic platform over a stationary ob- 


ject space field. Photograph 1 is taken at X q = X^ and yields 

11 2 

image points x , y . Photograph 2 is taken at X = X and yields 
P P o o 

2 2 

image points x , y for the same object space field X . The plat- 
P P P 

form is assumed to move with a three-dimensional velocity V. It 
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is known. from the basic lav/s of mechanics that X 2 - X" 5, = V«&t, 

o o 

2 2 

hence alternative expressions for x p , y p are: 

x 2 - x = -fflLCX - X 1 - VAt) / M, (X - X 1 - VAt) 

po Llpo ' 3 p o J 

y 2 - y n = -f [m„(x : - xj; - VAt) / - x* - v^t)] 

p o L 2 p o 3 p o -* 


(2j)2) 


The displacement of the image due to the movement between the 
two photostations may be computed as: 

A x ' x p - X P = - f k (x p - *0 - vat) / V X P - X o - T4t) ] 

- f k (x P - / M j (x P - 

•; C2.3) 

A y = yj' - yj = -f[M 2 ttp - X o * vat) / V x p " “ V4t)] 

- f [ M 2 (X p - X o> / M 3 (x p - fy] 

Expansion and simplification of equation set (2.3) yields the 

following image displacements: 

IL f (X - X 1 )*^ - VM_(X -X 1 )]^ 
liP o 3 3 P o_J 


Ax = -f^-2 i-E- _ 

M 3 (X p “ X o. V VAt).M 3 (X p - Xj) 


Ay 


M, fte - X^)M V - VM, (X - xh]l 
= _ f '-21 P 0 3__ 3 P ^_oJ_ 

M,(X - X 1 - VAt) *M (X - X 1 ) 

3 p o 3 p o 


( 2 . 4 ) 


Logically, image velocities due to translation of the photo- 
station may be readily computed from the above expressions as 


, ax 
x _ At 


v = &Z 

y At 


and 
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Similarly, it is possible to describe image changes due to 
changes in the attitude rates. It is recognized that 
\ 

M = M 2 ‘ VA 

M 3 

Then, by virtue of the double angle formulae, it is possible to 
define 


N = N_ . = ELR. ..R.R. , R R. . 

2 K *»At <p » «>*At 

W, 


Thus 


N + «»ht) R ( <f> + 4*4t) J5 ( w > + 


•where A = rate of change of k , 

$ = rate of change of £ , 

and = rate of change of <m . 


(See Appendix A) 


Hence the general expressions for image space coordinates of a 
moving system at any time with respect to a defined epoch, 

At = 0 are: 


x - x = -f [N. (X - X - VAt) / N_a - X - VAt)l 
po Llpo. 3po -* 


- y = -ffw_(x - X - VAt) / N,(X - X - VAt)] 
J o L 2 p o 3 p o J 


(2.5) 


In the above derivations, the rate functions have been assumed to 
be constant. However, this is a mathematical simplicity that can 
quickly be removed. 

By definition, image velocity is 5 and y and hence differ- 

3 ? 3 ? 

entiation of equation set (2.3) with respect to At yields the 
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following instantaneous image velocity expressions: 



y = -ff - X - VAt) - N_vl * [N,a - - V&t) 

P l L 2 p o 2j|,3p o J 


- [V X P - X o.- Vat) - V] '[ N 2 (X E - X o - V4t> ]}' 

• h (x P - x o - vat) J ' 2 

. dN- ■ 

where , similarly for ^ and N^. (See Appendix-'-A. ) 

Equation sets (2.5) and (2,6) are completely general and 
without restriction as to the initial orientation of the camera, 
the form of the object space, the time interval, At, over which 
computations are to be considered, 'or the magnitudes of the im- 
pressed rates. They may therefore be called "general" equations 
for the context of this study. 

2.2 Comparison of the General Theory with Kawachi* s Theory 

It can be demonstrated numerically that under the same ini- 
tial conditions similar values are obtained for the general theory 
and the pertinent formulae of Kawachi. However, a stronger case 

exists by virtue of the fact that It is possible to decompose 

15, 16 

equation set (2.6) into the forms given by Kawachi . This 

decomposition is demonstrated in Appendix £ using specific 
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examples. 

Thus, equation set (2.6) appears to be the general equation 
set for all forms and causes of image motion for frame-type 
cameras where the collinearity condition is the same as that 

initially assumed. The equation set (2.6) must be modified for 
panoramic photography, since equation set. (2.1) is not applicable. 

2.3 Some Implications of the General Theory 

Equation sets (2.5) and (2.6) were derived without any assump- 
tions as to camera type (other than those satisfying the basic 
equation set (2.1)) They are therefore applicable to between-the- 
lens shutter systems as well as focal plane shutter systems. In 
the former case, every image point moves during the time interval, 
At, that the shutter is open, whereas in the case of focal plane 
shutters each segment may be viewed as an independent exposure 
during which image blur may or may not occur. If the exposure of 
each strip is sufficiently short then no detectable image motion 
will be seen, although image distortion will be present. This 
distortion is fully explained by equation set (2.5) 

2.31 The 80 mm System of Lunar Orbit er 

The 80 mm system of Lunar Orbiter was a between-the-lens 
shutter system with one component of IMC driven from a V/H sensor 
designed to detect net forward velocity. Thus the film was ad- 
vanced or retarded at the required rate to minimize image blur. 
Variations in this rate have little meaning since it is a between- 
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the-lens system. The net effect of IMC is to produce a stationary 

photostation and consequently no photogrammetric reduction prob- 
lems are envisaged except those resulting from data transmission 
or from a lack of stability of the interior orientation elements. 

2.32 The 610 mm System of - -Lunar Orb iter 

The 610 mm system for Lunar Orbiter had a focal plane shutter 
adjustable for 3 exposure times ( 1 / 25 » 1 / 50 , and 1/100 sec.). 

IMC could only be- applied in the direction of flight and was com- 
puted from an on-board V/H sensor. Thus, the application of IMC 
is not constant unless the V/H ratio is constant. The photograph 
therefore suffers from two distortions: 

(a) That due to the focal plane shutter and predicted by 
equation set (2.5) if At of each segment from some epoch is known. 

(b) That due to forward motion IMC. This could be removed 
from the image coordinates on' a summation ttasis by multiplying 
the applied IMC rate by A the time lapse’ ^from epoch to moment 
of consideration. 

In the case of Lunar. Orbiters I, .II, III and V, both (a) and 
(b) type distortions are present, while for Lunar Orbiter I? ohly 
type (a) distortions are present.. (Actually,; some- blur; .is also 
present due to lack of IMC. This is becaus'e the direction”of 
appliable IMC was not in the direction of the net forward motion.) 
That is, for focal plane cameras in which IMC was not used, 
equation set (2.5) fully describes the relationships- between the 
object and image spaces. The important parameter is At, which 
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represents the segment exposure time along the x-axis with respect 

to some defined epoch, preferably defined as At = 0 when x =0. 

P 

A suitable method of obtaining At is available from the following 
data: 

(a) Exposure interval - controls blind slit width. 

(b) Blind velocity and direction. 

However as the blind slit width approaches the dimensions of the 
format, then the focal plane shutter system becomes equivalent 
to the between-the-lens shutter system. 

2.4 Some Experimental Results for Image Motion 

Using equation set (2.6) a number of computer runs were made 
simulating the 80 mm and 610 mm systems of the Lunar Orbiter pro- 
gram. The results of some of these simulations have been collated 
into Plates I and II; Plate I considers translational velocities 
and their effects while Plate II considers the rotational veloc- 
ities. both studies are made using the assumption that the photo- 
graphic system was initially vertical. Each plate has four major 
sub-sections dealing with each system at two different flying 
elevations. The sub-sections are composed of four separate graphs 
which indicate the expected magnitudes of the x- and y-directed 
image motion velocities for four important positions of the format. 
These positions are as follows, top to bottom and left to right, 
within each sub-section; 

(a) x = 0, y = max. 

P 

(b) x = max, y = max, an extreme corner. 

P P 
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(c) x = 0, v =0, center of format. 

P P * 

(d) x = max, y =0. 

P J p 

It is apparent from Plates I and II that not all image velocities 
are plotted. It may be assumed that those not plotted are zero 
or of such small magnitude that they can be considered as' zero. 

An analysis of Plates I and II readily shows that for rota-? 
tional velocities it is the focal length of the camera system 
that is important, regardless of the scale of the photograph. In 
the case of translational velocities, it is the scale of the photo- 
graph that is important in determining the magnitude of the ex- 
pected image veolcity. Unfortunately, long focal lengths are 
normally used to increase photographic scale and hence bath effects 
are present in the Lunar Orbiter systems, I-t is also most unlikely 
that the image motion effects can be made to cancel each other. 

2 .5 Conclusions 

Equation set (2.5) is therefore an adequate representation of 
the' processes controlling the functioning of both the 80 mm and 
610 mm systems of Lunar Orbiter IV. Should photography from the 
610 mm system of the other missions be considered for triangula- 
tion, then equation set (2,6) must also be considered. This 
would be quite complex, unless some simplifying assumptions could 
be made. 

It is suspected from previous experiences that the velocity 
terms are heavily correlated with their respective photostation 
parameters. If this correlation is of such a nature as to make 
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the system "singular", then it will be necessary to apply con- 
ditions on the velocity parameters in order to obtain acceptable 
solutions for the unknown parameters of equation set (2.5) « 

Equation set (2.6) fully describes the nature and magnitude 
of any expected image motion due to translation and rotation of 
the camera during the exposure interval. It therefore represents 
two functional relationships which describe the nature of image 
velocity in terms of the elements of interior and exterior 
orientation. 



3. MATRIX THEORY 


3.0 Introduction 

Methods for solving systems of normal equations date from the 
Gauss-Legendre era of the early nineteenth century, although the 
solution of small sets of simultaneous equations was already well- 
developed at that time. Thus, before the days of even the most 
modest digital computer, algorithms for either inverting or solv- 
ing small- to medium-sized systems were well-known. Since the ad- 
vent of the modern digital computer, the size of such systems has 
continued to increase, offering a unique field of endeavour for 
some scientists. 

Geodetic scientists, in their quest to accurately describe the 
earth and her near neighbours, have consistently been in the fore- 
front as users and developers of such systems. However, such use 
and development has seldom been tempered with an adequate numeri- 
cal analysis of the problem. 

In photogrammetry , the problem of inverting and/or solving a 
large system of equations resulting from a simultaneous, multi- 
station triangulation is already formidable. The proposed augmen- 
tation of the collinearity conditions places an even greater bur- 
den on the algorithms and computing machinery in present use. 
Furthermore, it has long been recognised that it is not normally 
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possible to simultaneously isolate both the camera constant, f, 
and the elevation of the exposure station, Z Q , as these para- 
meters are strongly correlated. This lack of separation can be 
shown by investigating the ,, condition ,, of the matrix. 

In this chapter, therefore, methods for detecting and over- 
coming the problems of size and stability will be discussed. 

3*01 Review of Matrix Inversion 

18 

Most texts on matrix algebra, e.g. Faddeev and Faddeeva , 
not only provide a rigorous explanation of the matrix inverse, but 
also provide numerous algorithms for computing the same. Further- 
more, program libraries such as the IBM Scientific Subroutine 
19 

Package already provide the photogrammetrist with the necessary 
computer software for many of these algorithms, thereby greatly 
reducing his work. 

However, once the core storage capacity of the computer is 

exceeded, then the photogrammetrist must again address himself to 

the problem of devising adequate computer software. Consider the 
T —1 

matrix system (A. £ A) of normal equations. Classically, this 

is a banded system. However, in the event of either A or 2 

T —1 —1 

becoming full, then A Z A is full. The case where ."5J is full 

rather than block diagonal is especially realistic when satellite 

orbital data are used to constrain the photostation coordinates. 
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Table 2 lists some important references pertaining to the 
solution of (A 2T A,) - . Figure 2 depicts diagrammaticaliy the 
form of the normal equation matrix. 


Table 2 

Reference List Pertaining to the Solution of High Order Inverses 


Matrix 

Type 

Method Name 

„ 

Reference 

1 

Partitioned Regression* 

Brown 20 ’ 21 " 

2 

Triple Block Method* 

22 

Snowden 

2 

Elassal’s General Algorithm* 

Elassal 25 

3 

Triangularization 

Uotila 24 , IBM 19 

4 

Gauss- Jordan 

Berezin and 



Zhidkov 29 , 
IBM 19 

5 

Partitioning 

Faddeev and 


- 

Faddeeva 1 ^ 

6 

Successive Partitioning 

Snowden 22 , 
Berezin and 
Zhidkov 29 


"■These methods use the same fundamental concepts and 
essentially differ only in name. 









figure 2 

Some Common Normal Equation Forms 
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m 

. xxxxxxxx 

XXX 

Banded Matrix, 

- 

XX 

XXX 

Banded Matrix, 

, xxxx 

often solved by 

XXX 

often solved by 

X XXX 

Partitioned 

XXX 

Partitioned 

X XXX 

Regression. 

XXX 

Regression, 

X XXX 
X XXX 

X XX 

m « 


XXX 

XXX 

XX 

• « 

Triple Block, 
or Elassal's 
Algorithm. 

“xxxxxxxx' 


“ xxxxxxxx ' 


xxxxxxx 

Symmetric 

xxxxxxxx 

Full Matrix, 

xxxxxx 

Matrix, often 

xxxxxxxx 

often solved by 

xxxxx 

solved by Tri- 

xxxxxxxx 

Triangulariza- 

xxxx 

angularization , 

xxxxxxxx 

tion, or Gauss- 

XXX 

or Gauss-Jordan. 

xxxxxxxx 

Jordan. 

XX 

-X 

. 


xxxxxxxx 

xxxxxxxx 

# 


* xxxfxxxjxx" 


*. ■ 

xxxxxxj; x 


xxxjxxxxx 

Partitioned 

xxxxxxx 

Positive Semi- 

xxxixxx-xx 

Full Matrix 

“xxxxxx 

Definite Matrix, 

Txxxxxjxx 

illustrating 

xxxxx 

often solved by 

xxxxxxxx 

Successive 

xxxx 

Triangulariza- 

xxxxxxxx 

Partitioning 

XXX 

tion, or Gauss- 

XXXXXXXX 

Scheme. 

XX 

Jordan. 

XXXXXXXX 

« I 


L 2 . 



It is recognized that Table 2 is not complete and that the 

literature contains many other algorithms. One such algorithm 

that appears to be growing in importance is that used to "border" 

a matrix. This technique allows the use of the positive definite 

character of the matrix to the stage where the matrix becomes 

26 

positive semi-definite (see Needham ) . 
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J.02 Review of Solution Methods 

As in the case of matrix inversion, a wealth of information 

can already be found in standard texts, e.g. Faddeev and 
•18 

Faddeeva , and in the journals. Unfortunately, most methods of 
solution depend on the normal equation matrix being positive 
definite, a feature that is destroyed when constraints are added 
to the system. However, these methods have usually proved to be 
superior to inverse techniques in both time and stability. 
Unfortunately, the variance-covariance matrix of the adjusted 
parameters is not easily obtained. 

A complete description of the following methods can be found 

18 _ _ . T5 -* j rr 1 25 . 


in either Faddeev and Faddeeva 
Gauss Elimination 
Square Root Method 
Gauss-Seidel . 

Relaxation 


1 


or Berezin and Zhidkov 

These methods are normally used on 
symmetric- and full-type matrices, 
rather, than oh banded systems. 


3-1 Vector and Matrix Forms 

Since the concept of a vector norm is more easily envisaged 
and understood than that of a matrix norm, it is considered 
logical to develop the background for the n x 1 vector X before 
considering the more general matrix A of size n x m. 

The norm of a vector is defined as a real, non-negative 
number, expressed as Jjxjj which represents, in some manner, the 
size of the vector X. The following three manipulative rules are 
important : 
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It xj| *^= 0 ( || 0 j = 0, otherwise ;?• 0) 

||kx|| = ||k|«|x| (3.1) 

|)x + Yj « Jx|| + Jy|J 

Consider the following special cases for the n x 1 vector X: 

(a) n = 1. In this case the size of X is best expressed by 
the modulus jxj . However} it is not the only estimator. 

(b) n = 2. There is now no single number that accurately 
gives the size of such a vector. Some of the most logical and 
significant estimators are: 


C I a K| 


+ X 



(c) n = 3* As the size of the vector increases, so does the 
number of logical and significant estimators. Some of these es- 
timators are as follows: 





From these examples it is readily seen that a general defini- 
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tion, the Holder norm , is possible* 

Ixl 


X 


X . 
X 


1/k 


(3-2) 
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The term Euclidean norm is applied to the k = 2 norm, since 
it represents the length in multidimensional space from the point 
to the coordinate origin. 

From the Holder norm and from .the above examples, it may 
readily be deduced that the following inequalities exist between 
the norms: 


*1* 


*II X IU 


2 


- X 




CO 


It is now considered pertinent to consider norms associated 
with matrices and systems of linear equations. The basic alge- 
braic rules which follow are similar to those for vector norms 
(equation set 3.1): 


A 


kA 


3* 0 ( ||o 


= 0, otherwise 0) 


= k 


A + B « 


AB 


B 


( 3 . 3 ) 


B 


AX 


18 


Faddeev and Faddeeva have shown that the following forms corres- 
pond to the equivalent vector norms. Thus, if the second norm is 
being used in vector work, then the associated second matrix norm 

must be used for associated matrix work. 

n 


= max 


■<E 

i=l 


a . . 
13 


) - a column norm 


T T 

= max ()0 , V of the matrix A A - A must be Hermitian 


= :'.ax ( 2 
i j=l 


a. . ) - a row norm 

10 
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It is obvious that the computational labour involved in 
the spectral norm, is immense. Therefore it has been common to 
use the more easily computed Euclidean norm, which is generally 
greater than the spectral norm by a factor of up to ifn. 

The Euclidean norm is 

34 (3-4) 


2 ’ 



1 

r n n . 

“ 

A l 

E 25 

£ £!■>. 

2 


The following identities hold for matrix norms : 



Two very common matrix norms are defined as follows; 

(3-5) 


M(A) j= n max 

1,0 


a. . 
ID 


N(A) = 


3-D 


a. . 
3-D 


% 


= trace A^A (3.6) 


It is evident' that the norm M(A) belongs to the k = co class, 
while N(A) belongs to the k = 2 class. Some important relation- 
ships between norms are now given; 


J M(A) $ 


“ M(A) 4 || A"j 


~ M(A) 4 


< M(A) 
L <? M(A) 
, 4 M(A) 


M(A) 4 K(A) 4 M(A) 


i N (A)*||A 
| »<A) < 1 A 


2 4 N(A) 

4 Vn N(A) 
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^?:U) £ || a|| x ^ k(a) 

sHIa # | A ||» 4 ■F Ha 

H4 s || a ||i * & Ha 

kII a MHi «"H- 

1 (I + C) 1 4 m 

I 1 - (I + crl || < iqff 

The last two inequalities are proved in Appendix C , 
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orthogonal matrices in which at least P = H = 1. Unfortunately, 
the upper limit is not bounded and therefore experience in 
working with systems becomes an important part of the decision 
process. Much of this experience is gained by experimenting with 
test matrices. 

3*21 Orthogonal Test Matrices 

It is possible to construct orthogonal matrices of any de- 
sired size for software testing. It is recognized that inverses 
computed via algorithms which ignore orthogonality should ex- 
perience minimal roundoff error. This is due to the low magni- 
tude of the stability numbers. However, such matrices seldom 

.occur in Geodetic Science. A useful test form is given by 
30 

Newman and Todd : 

A = ( a . . ) where a . . = , — ~ * sin f - ) 

v 13 ' aj Vn+1 v n+1 ' 

and since A is orthogonal 



3*22 Ill-Conditioned Test Matrices 

Just as it was possible to construct matrices with highly 

stable characteristics, it is also possible to construct matrices 

with highly unstable characteristics. A finite segment of the 

well-known Hilbert matrix and its variants have been shown to be 
30 31 32 

highly unstable"^ ’ y J . For instance, the P number associated 

3*5 n 

with these matrices is approximately e , where n is the order 
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of the matrix. Thus, this class of matrices, while not repre- 
sentative of the type of matrix normally encountered, offers the 
possibility of comparison between two given algorithms. The 
form of this matrix is as follows: 


then A 


= where 

-1 


a . . = 


i+j-1 


= J 


(-l) 1+3 (n+i-l)i (n+j-1)! 


where b. . s= 


(i+j-1) [(i-l)i C 2 (n-i)J (n-j)l 


3.23 Some Tridiagonal Forms 

Tridiagonal forms of the normal equation matrix do exist in 
Geodetic Science ^ . (For an example see Appendix D.) Further- 
more, many tridiagonal forms have known eigenvalues and hence 
known P numbers, in addition to an algebraic inverse. It is this 
class of matrices which should be used to investigate computer 
software performance for geodetic applications. Gregory and 
Harney list a number of suitable forms. However, experience 
has shown that the computational merits of a particular algorithm 
are not evident up to an order of at least n = 75 • The reason for 
this is that the test matrices were positive definite and 
possessed uniformly small elements which could be exactly repre- 
sented. Under such conditions, many algorithms will yield 
similar results, 

e.g. for n = 70 and a tridiagonal Matrix given by 

A = fa. ,\ where a. . = 2, i = j 
1 xj; xj j 

&i j = 1 . I X- j | = 1 

\ 3 = °’ l i - U >. 1 



35 


then A -1 = (b^ = S+I^ c i^ where c i; j = i(n-i+l), i = j 


c . . = c . . -i, i > i 

i,a-i ’ J 

C ij “ C ji, j < i 


Four representative methods yielded the following values for 


Desk calculator 
DSINV (IBM 1 9) 
DMIl'IV 


11 *' 

0.985 915 492 957 746 5 
0.985 915 492 957 748 
0.985 915 492 957 749 
VERSOL (See Appendix. E) O .985 915 492 957 748 
The P number for this problem is approximately 1985 * Thus- 
it is seen that,'- for this class and size of matrix, these inverse 
subroutines yield results with minimal roundoff error, despite 
the large P number. 


9 


3*24 A Test Matrix 

A large matrix (370 x 370) resulting from the adjustment of 
gravity observations was available for testing. This matrix had 
gravity differences for its basic data. Thus, while the matrix 
was formed without any regard to a pattern, it could easily be 
condensed into a variant of the tridiagonal matrix. This matrix 
therefore became a convenient test matrix, as it was- assumed to 
be a more representative matrix than' the theoretical tridiagonal 
case due to the size of its elements as well as its order. 

Since the matrix has an impressed condition, it was capable 
of being reduced -into two variants.- Table 3 .shows vfthe -stability 
indicators for both variants. As mentioned earlier, ho definitive 
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statements can be made regarding the quality of the matrix from 
these indicators, since they possess no upper bound. However, 
experience indicates that an acceptable solution to the inverse 
problem was obtained. 


Table 3 

Stability Indicators for the Test Matrix 


Indicator 

Variant 1 
n = 369 

(no imposed condition) 

Variant Z 

n = 370 

(condition imposed) 

M 

^3.1 x 10 

^ 9.2 x 10 6 

N 

?“ 3.4 x 10 3 

> 1.2 x 10 4 

P 

4.2 x IQ 4 

5*5 x 10^ 

' 


3.3 Perturbation Theory 

It was shown in Section 3*2 that the classical indicators are 

of limited use since they have no supremum. For this reason, 

other indicators of stability have been sought. A successful 

method based on the exact computation of a perturbed system is 
33 

due to Wilkinson . The method seeks to place bounds on the per- 
turbations of the system necessary to perform a computation, rather 
than to follow the forward error propagation. 

Consider the following matrix system: 

AX + Y = 0 (3.H) 

assuming that A is square with a non-zero determinant. 
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Then if A is perturbed by an amount E, X is perturbed by the 

vector,$X such that the equation of the system becomes 

(A + E) (X- + &X) = -Y = Z (3*12) 

or A(I + C) (X + Sx) = Z (3-13) 

-1 -1 

where AC = E, which implies that, since A exists, C = A E. 
Alternatively, |'cj 1 is a sufficient condition (See Appendix C) 
Solving for &X yields 


Sx = (I + C) _1 A -1 Z - X 


,-lr 


but from equation (3*11) X= A Z 
Hence Sx = [(I + C)” 1 - i] X 

Application of the norm theory to this expression yields 

HI - ||[ (I + c)_1 - jJ*! 16 ||[ (I ♦ c) " 1 - I ]|-J?l 

Sx|| 4, [|(i . O)- 1 ! , [-ij]-||x| 

8§x 


u i 


. 1 + (1 + c) 


,-i 


I 

But C = A“ 1 E 
Hence IjSxll 

"W 


Hell 


A -1 E 


l-llA^E 




,-1 


E 


( 3 . 1*0 


1- A 


,*•1! 


E 


This last equation expresses, in terms of norms, the relative 
change in a solution vector due to a perturbation in the original 
matrix A, of size E. It must be stressed that the above expres- 
sion is an inequality and that only an upper bound has been deter- 
mined. This bound could be in considerable error if 
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It is important to note that the condition, number is the 


decisive quantity in determining an upper bound. Consider a con- 


dition number defined as 


k = lA • 


.'-1 


(note M = -•H(A)*M'(A _1 ) and N = i N(A) -WCa" 1 ) ) 
n n 

and a relative perturbation S defined as S = jrff- 


Then the relative change can be expressed as 


JCIL 

«*i 


f — t 

J 

< 

• a| • 

|e! 

/ 

[l-- I A- 1 

• 

A l‘ 

bi 


II 

w_ 

/ 



11 ‘ 

n 


fi 


= kS 


1-kS 


(3.15) 


Ordinarily $ Is small, but if kS-yl due to large k, 
then |Sxil -^-oo • This is critical in digital computers where 


Sis a fixed ratio. 

2 ? 

It is possible to extend this perturbation theory to include 
perturbation on the constant vector Y . Thus the- perturbation 
equation (3.12) becomes 


(A + E) (X +$X) = (Z«+ SZ) (3.16) 
Application of the above theory yields the following equivalent 
expression to equation (3.1^). Note | $»zj| = |Sy| 



(3.17) 


The most efficient method of estimating |Sx| from a computa- 
tional viewpoint is with the co norm. The sacrifices in accuracy 
made with this norm can be determined, if desired, from the 
identities of Section 3«1» However, the physical significance of 
this norm is more important. Thus equation (3«17) becomes 
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Sx 


Co 




A _1 oo * 

e IL* Ik L + 

Ik 

1 L-1syL 

1 - | 

A’ 1 

• 

oo 

E 

f 

CO 


(3.18) 


It is to be noted that the following definitions apply to 
infinity norms which can be physically interpreted as row norms, 
since the maximum element vector norm is a degenerate case of the 
more general matrix row norm; 


lx 


= max x. 
co x| 


n 


= max 


co 


(2 k, ) 


3=1 


1 3 


- maximum element norm 


- row norm 


The perturbation introduced in a digital computer when storing a 

number is given by 1:2 where b is the number of bits in a word. 

For the IBM 360/75 of The Ohio State University Instruction and 

Research Computer Center (IRCC) the perturbation level is 

7 

approximately 1 part in 10 for single precision and 1 part in 
lO 1 ^ for double precision. Thus, |e|^ aDd |SYj|^ are governed by 
the size of the elements of A and Y. For this discussion, it is 
convenient to view the perturbations as decimals. This is the 
situation which would exist if row and column normalization of the 
system were performed prior to solution. It is also convenient to 
use the double precision mode figures. Hence 


EIL = n-10 


-15 


|sy| = 

11 *ioo 


10 


-15 


For the test matrix (Section 3.24), the following values of 


A|| and X were obtained: 

£* CO 



bo 

» = 369, jA- 1 ^ = 1.2 

n = 370, ||a- 1 L = 44 

ML - 200 

When these values are substituted into equation (3«l8) together 
with the values of | E and Yj] ^ the following are obtained: 

n = 369 , | SxJ ^ 10“^° 
n = 370, |$X| c- 10“ 8 

‘This indicates that the solution vector is accurate to 10 
significant figures for the n = 369 configuration and 8 signifi- 
cant figures for the n = 370 configuration. It therefore appears 
that the added condition did not strengthen the stability of the 
matrix. Moreover, the stability of the matrix is not a direct 
function of the P, N, and M numbers when the matrix is of a high 
order . 


3.4 Matrix Refinement 

This method, which is based on the control computation AA \ 
seeks to take an inverse which is known to possess roundoff 


errors and refine it until these errors are no longer of any com- 
putational consequence. The discrepancy between AA and I is 
naturally an indicator of the degree of stability, for the closer 
AA -1 is to I, the more stable is the computation. Consider the 


following definition of the discrepancy matrix: 



It is immediately recognized that for A Q 
approximation to A \ then jc^l £ k <• 1. 


(3.-19) 

to be an acceptable 
The more easily com- 
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puted 1st and 00 order norms are most commonly used, although 

Hotelling‘S introduced the concept of refinement using N(A). 

3*3 

Faddeev and Faddeeva show that the following relationship 
exists : 


A*" 1 = A’ 1 , (I + C , ) 
m m-1 m-1 


Hence, from equation (3.19)t 

C = I - AA*" 1 
m m 

Expanding yields 

C= I - AA” 1 1 (I + C ) 
m m-± m-1 


= i - (i - c Jd + c j 

m-1 m— 1 


( 3 . 20 ) 


- 2 4 

= cr _ = c _ 

m-1 m-2 




Hence, 


m 

A * = A“ X (I - CZ ) 
m O 


(3.21) 


-,m 


2 “ -1 -1 

and in the limit, as — 3»- 0, A‘ — s- A 

0 m 


This may be re-arranged to yield the following equation on which 
the computational algorithm has been built : 

( 3 . 22 ) 


-1 -1 -1 -1 

A = 2A . - A \AA - 
m m-1 m-1 m-1 


It is now desired 'that a single number, rather than the dis- 
crepancy matrix C, should indicate the quality of - the computations. 
Thus, the concept of the norm must again be introduced. 

Consider 

D = A’ 1 - A -1 
m m 




r- 

1 

r*i 

1 



!! 

A - A 

m 

ti 

- A C 0 


then 
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(3.23) 


( 3 * 24 ) 


The test matrix of Section 3*24 was subjected to this refine- 
ment process. The initial and final values associated with the 
test matrix are listed in Table 4. The matrix was subjected to 
two refinement steps. 


Table 4 1 • 

Some Values for the Test Matrix Associated with Refinement 




C 0 


m ■ 

1 



Initial 

Final 

Initial 

Final 

Variant 1 

n = 369 

4 x 10 -11 

1 x 10“ 10 

„ --20 
2 x 10 

2 x 10~ 20 

Variant 2 
n = 370 

1.2 x 10~ 3 

2 x 10" 10 

1.4 x 10~ 2 

if x 10~ 18 


It is evident from Table 4 that only Variant 2 benefited 
from the refinement process. That is, the refinement process 
roundoff errors were themselves the limiting factor in the 
369 x 369 case. However the 370 x 370 matrix, which was inverted 
by using the bordering technique on row 370, did respond to re- 
finement. This response was due to the fact that the bordering 
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algorithm produced roundoff errors in excess of those restricting 


the refinement process. 

A number of test examples of small to moderate order were 
tested, including finite segments of the Hilbert matrix. Results 

from .these tests indicate that refinement is not possible for P 

5 

levels in excess of 10 , and that .system instability occurs for P 
levels greater than 10^. 
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These findings are consistent with those of Wilkinson , who 
shows that inequalities of the following type exist. (Note: With 
these inequalities, the quality of the initial approximation must 
be carefully considered.) 

(a) Refinement occurs provided 

I A" 1 2 k-1 /n < 2 -1 

i oo 

where k is the number of machine bits in a word. 

(b) The gain in precision in binary bits per iteration, m, is 

n2 ' k |ML< 2 ' m 

(c) Component error level for residuals is approximately 



where 2 ? is the maximum component of A 

(Note: For first iteration, the level isv^ times above.) 

(d) Error level in solution given by 

(n)* HI ®] 1 

where i is the iteration. 

Three examples illustrating these computations are listed in 
Table 5« These computations are pertinent to the IBM 360/75 of 
IRCC, where k = 5&, hence ^ 7.6 x 10^ and 2 ^ ~ 1.3 x 10 ^ . 



Table 5 


Inequality Tests to Determine Possibility of Refinement 

























Table 5 indicates that any instability in the. systems would 
not become apparent until the second iteration if the initial 
inverse choice were good, since the errors present are not of suf- 
ficient size to unduly perturb the inverse. 

However, when these errors are augmented with the larger 
errors of the refinement process, then an inverse quickly ceases 
to exist. When the. refinement errors are of moderate size, or 
comparable to the roundoff error, oscillation takes place. Such 
an oscillating system was observed for the n = 6 segment of the 
Hilbert matrix. 

Unfortunately, this technique may be impractical except in 
special circumstances, as central processing unit (CPU) time was 
33 minutes per iteration for the 370 x 370 matrix. This time 
would certainly decrease with higher rates of information transfer 
between the disk unit and the core. The rate of transfer used in 
this problem was 312K bytes per second. Additionally, if a greater 
core region was used, then fewer calls to the IBCOM routines 
would also reduce CPU time. The present program required 16K 
bytes of core for storing the associated instructions. The neces- 
sary four work vectors require additional space which in this case 
amounted to l6K bytes. This space is determined by the order of 
the matrix. 

The method could be useful in extending the core range of the 
machine, since it is theoretically possible to obtain a single 
precision inverse of limited quality, and thence to refine it to 
the desired level. However, care must be exercised to ensure that 
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the inverse is acceptable. 

The program for accomplishing this work is listed in Appendix 

E. 


3*5 Reinforcement 

The principle of reinforcement is that the matrix A is con- 
sidered as the last term of the following sequence:. 1 


A 0 " I * A l* A 2 *** ^k-l* \ 


A = A 
n 


That is, matrix A^ is obtained by replacing the kth row of A^._^ 
with the kth row of A^,, thereby building up the desired matrix 
which at this juncture is unspecified. 

zl3 

The following derivation is due to Faddeev and Faddeeva : 
Consider the matrix A to be non-singular and with known inverse, 
and consider the column vectors U and V defined such that 


UV = 


* 


• 


m 

U 1 


U 1 V 1 

U 1 V 2 

... u^v n 

U 2 

• 

• 

[V v 2 "• V J = 

U 2 V 1 

U 2 V 2 

... u_v 

2 n 

♦ 

u 


U V- 

u v_ 

. . . U V 

n 


n 1 

n 2 

n n 



» 


. 


(3*25) 


Clearly, UV has rank 1, since every row of equation (3*25) is a 
linear combination of another row. 
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Then it has been shown by Dwyer and Waugh that for the 


matrix 

B = A + UV 

B" 1 = A" 1 A -1 UVA _1 

1+VA U 

provided that 1 + VA ^U ^ 0 


(3*26) 

(3*2?) 



hi 

That is, it is possible to find the inverse of B, given a matrix 


A with known inverse which differs from B by a matrix of rank 1, 
In particular, if the UV matrix is constructed as 


= ° l V l v 2 •” V k ••• %] 


then only the elements of the kth row are being changed. 
Using these concepts, equation (3.2?) becomes 


B - 1 = A ’ 1 


— — a, (VA*" 1 ) 
! + Va k 


( 3 . 28 ) 


where a, is the kth column of A . 
k 

It is now recognised that VA is a row vector resulting from a 
summation over a column. Hence, the following holds: 


a fc (Va ) 


b . = a . 

3 3 (l+Va ) 


(3-29) 


-1 -1 

where b . and a . are the jth columns of A and B respectively. 
I 0 

Finally, the series concept is added, which provides a known 
inverse for A. 


Thus, equation (3.29) becomes 

, (k-1) _ (k-1) 
„ b (k-l) . \ Yb .i 

3 3 (k-1) 

l+Vb. 1 * 1} 
k 


(3.30) 


v/hich is the working algorithm for the solution of the inverse. 

(k-l) 

It should be -noted that the only B for which an inverse 
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is known is the identity matrix I . This matrix differs from the 

given matrix. A, by n rows of the form (v v . . . v, ... v ). 

jl c. n 

Hence, the kth vector of v's is the kth row of the given matrix, 
with the exception that the kth element must be zero to maintain 
the desired linear relationships. Thus, 

V becomes a., ... a^.O, a fc+1 , ... 

It therefore becomes necessary to pre-multiply A by E2~type 
elementary transformations to achieve unit diagonal terms to 
satisfy the above conditions. This, in turn, necessitates post- 
multiplication of B by E^-tyP® transformations to return the de- 
sired inverse. It is also noted that the method applies specifi- 
cally to positive definite forms, which is a considerable draw- 
back. 

A computer program utilizing the above method is given in 
Appendix E. This program makes use of the disk storage and there- 
fore may -be used to invert any full positive definite matrix, 
since only six vectors of the order of the matrix are needed to 
accomplish the inversion. 

Unfortunately, the method is very slow. Inversion of the 
369 x 369 test matrix required 210 minutes of CPU time. Moreover, 
the number of significant digits obtained, when compared to the 
refined solution, was such that the process could be termed "un- 
stable". Only four significant digits were obtained. Thus, it 
seems that until faster transfer rates can be realized together 
with increased word bit size, this method must be passed over as 
a mathematical curiosity. 
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3.6 Solution by the Square Root Method 

This standard method of obtaining a solution vector is very 
stable for positive definite matrices and quite fast, computation- 
ally, on digital computers. Furthermore, since it is essentially 
a row process, it is readily adaptable for unlimited size by inter 
locking the auxiliary disk storage facility with the magnetic core 
It can be extended to positive semi-definite matrices by simply 
defining A as a Bermitian matrix, rather than as a real matrix. 
Furthermore, this change of definition need not be invoked compu- 
tationally until the positive definite character of. the matrix is 
destroyed. 

Consider the following matrix system: 


AX + Y = 0 


T 

Then A can be' triangularized such that A = BB 
where B is triangular. 

Then solve B*Z- = -Y by substitution for Z 
T 

and finally B X = Z by substitution to determine X, the 
solution vector. 

In the above. 


i-1 


B . (b > . 


..-Si 

T 1 "• *1 


11 


A 10 


* 


(i > 1) 


i-1 




kTl lk 


b. . 

11 



(j > i) 


0 


(i > 
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The main difficulty in this efficient procedure is that the 
variance-covariance matrix of the adjusted parameters is not de- 
termined. Thus, other methods of determining the inverse must be 
undertaken. 


5*7 The Monte Carlo Method 


... The Monte Carlo method may briefly be described as 
the devic'e of studying an artificial stochastic model of 
a physical or mathematical process. The device is cer- 
tainly not new. Moreover, the theory of stochastic proc- 
esses has been a subject of study for quite some time, 
and t\ie novelty lies rather in the suggestion that where 
an equation arising in a nonprobabilistic context demands 
a numerical solution not easily obtainable by standard 
numerical methods, there may exist a stochastic process 
with distributions cr parameters which satisfy the equa- 
tion, and it may actually be more efficient to construct 
a process and compute the statistics than to attempt to 
use those standard methods ... -q 

Householder' 5 


The Monte Carlo method appears to offer a convenient method 


of overcoming the lack of an inverse when a solution method is 
used. The theory and development of Oswald^ was converted for 
application to the IECC system, but as yet suitable and consistent 


results have not been obtained on medium scale test matrices. The 


reasons for this are likely to be many, but the problem appears to 
lie in the random number generator or its application. 

At present, several thousand walks are required for 3 and k 
digit' accuracies, thus giving rise to fairly lengthy execution 
times. A reduction in calls to the random number generator and 
better random numbers would materially aid in the solution of 
these problems. However, for the moment, this task does not 
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warrant further pursual, although it should not be forgotten. 

An interesting aspect of the Monte Carlo method is that if 
the finite variance condition is fulfilled, as it must be for an 
inverse to be obtained, then roundoff error plays an extremely 
small part in the computational process. Rather, the limiting 
factor appears to be the number of walks necessary to achieve a 
desired level for the elements of the 'variance-covariance matrix. 

In many cases, three significant digits would suffice. 

3.8 Conclusions 

In summary, it is to be noted that the inversion methods in- 
vestigated in this chapter do not appear to approach the speed or 
accuracy of those commonly used and mentioned in the review. How- 
ever, this should not foreclose the possibility of further work on 
this subject. In particular, methods for quickly and accurately 
obtaining elements of the variance-covariance matrix associated 
with a solution vector should receive attention. The Monte Carlo 
method is but one suitable method which may be able to yield 
suitable results in moderate CPU times. 

With regard to the problem of stability, it has been forcibly 
pointed out that the classical indicators are of limited use and, 
in general, it is far better to compute the norm of the solution 
vector and to compare this norm with the required precision levels. 
Should the precision of the solution be unacceptable, then refine- 
ment may be attempted, if feasible, or another solution tried. 
Alternatively, the computational stability of the matrix may be - 
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improved by applying suitable constraints to the normal equation 
set. This may upset the positive definite nature of the matrix, 
but this is of no great importance. Should the stability decrease, 
then critical consideration should be given to the desirability 
and necessity for the constraint, and to the efficiency of the 
computational process used. 



k . NUMERICAL TESTING 


4 .0 Introduction 

In Chapter 2 the augmented collinearity condition equations 
(2.5) were proposed. Chapter 3 developed a number of numerical 
concepts and tests to provide a rigorous mathematical basis for 
determining whether or not the collinearity conditions can be satis- 
factorily augmented with additional parameters. It is therefore 
the purpose of this chapter to link the previous two sections by 
means of numerical tests as well as by a clear demonstration of 
the advantages of the augmented system. 

Two principal methods of testing are available to research 
workers. The first method involves testing with simulated data. 

In essence, this method consists of mathematically generating ar- 

* 4 . . 

tificial point data according to some defined relationships, then 
using this data to test and check the new functions and associated 
computer software. The resulting estimated parameters may be com- 
pared against the true or known values of the parameters. 

The second method uses real data derived from the observation 
process. The parameters estimated in this way cannot be compared 
against known standard values, since these do not exist. The ul- 
timate object of the method is the determination of these unknown 
parameters to within a given confidence interval. The smaller the 
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given confidence interval, the more accurately is the unknown 
parameter determined. 

In general, the first method allows a greater variety of tests 
to be made. However, there is always a small probability that the 
simulated data conform to the proposed model, whereas the real 
data are not represented by the proposed model. Thus model testing 
is not complete without real data tests, even though tremendous 
insight and understanding may be gained when simulated data are 
used. Unfortunately, real data are often most difficult to obtain, 
and even more difficult to implement initially. The Lunar Orbiter 
IV mission is not an exception in these matters, as contractual 
difficulties have delayed the furnishing of such real data. Thus 
only simulated tests can be described in this report. 

The tests to be described in this chapter fall into two cate- 
gories. The first is concerned with single photo tests, while the 
second section builds the single photo tests into a multiple photo 
block. 

k .1 Generating the Simulated Data 

The well-known collinearity condition equations (2.1) have long 
been recognised as the mathematical expressions by which object 
space points are connected to image space points by way of the lens 
nodal points. It was therefore natural to choose these equations 
as the basis on which the simulated data were generated. 

The object space field was selected from the Department of 
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Defense Selenodetic Control System 1966^. It consisted of all 
points between 0°N and 30°N and between 10°W and 10°E. These 
points were then transformed into a right-handed X, Y, Z cartesian 
system. The selected object space field roughly corresponds to 
the area recorded on photographs numbered 109 and 102 of Lunar 0r- 
biter IV. Figure 3 illustrates the area- chosen and the limits 
of the photographic coverage. 

Photostation positions approximating the position of exposures 
109 and 102 were computed in the X, Y, Z system, together with 
attitudes in the omega, phi, kappa rotation system. It was there- 
fore possible to compute the corresponding image coordinates for 
each object point at each of the two assumed photostations for a 
610 mm focal length system. Figure 4 illustrates th* density 
and location of image points for a photograph similar to photo 
number 109. The same object space field was transformed into an 
image space for slightly different photostations corresponding to 
pseudo-velocities. Next, the image space was split into eleven 
equal sections corresponding to ant time interval of 0.01 sec. 

It was assumed that the focal plane shutter progressed from posi- 
tive y to negative y at a rate of 1 m/sec. The positions of those 
points which fell in a particular at interval were therefore com- 
puted for the photostation corresponding to that pseudo-velocity 
times At. The initial epoch. At = 0, was assumed to occur when 
the shutter crossed the midpoint of the image space. This is il- 
lustrated in Figure 4. 
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FIGURE 3 

LUNAR OttlTIt IV PHOTO INDEX 

(near Slot) 

WITH OUTlINf Of ARfA USED IN 
GENERATING THE SIMULATED DATA 




FIGURE 4 


TYPICAL SET OF SIMULATED DATA POINTS 
ILLUSTRATING THE UNIFORM DISTRIBUTION 
AND THE ELEVEN 0.01 SEC. EPOCHS INTO 
WHICH THE PHOTOGRAPH WAS DIVIDED 



X 


DIRECTION 
OF BLIND 
TRAVEL 


VJ1 

-o 
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Thus an image space field with a corresponding object space 
field was constructed without resorting to the proposed model, 
although the concepts are similar. As is the case with such simu- 
lated data, the determined parameters possessed true values against 
which they could be tested. The total number of points in the ob- 
ject test area was 46, 32 of which are shown in Figure 4. The 
exact number and location of points used varied from test to test, 
depending on the particular requirements of the test. 

4.2 Single Photo Tests - The Resection Problem 

Usually space resection problems are concerned only with the 
elements of exterior orientation. However, in the case of the 
Lunar Orbiter missions where the coordinates of the principal point 
must be considered as unknown, it is necessary to include these 
quantities amongst the unknown parameters or quantities to be ad- 
justed. Furthermore, there is evidence to suggest that these quan- 
tities are not fixed, unlike the case of frame cameras used in 
conjunction with recoverable film. It is therefore convenient for 
the purpose of this study to consider the unknown coordinates of 
the principal point as part of the exterior orientation elements. 

There are three principal subdivisions to the space resection 
problem. These subdivisions are made according to the amount of 
observational data that are available. They are as follows :- 

Case (a), characterized by - observations on photopoints. 
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Case (b), characterized, by 


r 


observations on photopoints, 
observations on elements of 
exterior orientation. 


Case (c), characterized by - observations on photopoints. 

- observations on elements of 
exterior orientation. 

- observations on survey points. 


4.21 Case (a) 

In accordance with the above description of this' case, the 
general mathematical system F(L ,X ) = 0 was chosen. This expres- 
sion is linearized according to the usual Taylor series expansion 
method to yield: 

AV + BA + £ =0 (4.1) 


In this case, 

<>F 


^-P 


cs x. 


A = 


<)F 


2£ 


^F 






^F 


2P 


)y r 


o 


Hence, (4.1) becomes 
E E 

V + BA+e = 0 


0 


= I 


1 


The system is completely general for n points. Under such con- 
ditions, the dimensions are as follows, 
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E E 

V + B • A + € = 0 (4.2) 

(2n x 1) (2n x 6) (6 x 1) (2n x 1) 

and the solution or correction vector is the well-known form 


E 

A = - 


T 


E 
P B 



e t 

B P 


e 

E 


(4.3) 


The correction vector A is now added to the approximate value 

vector X qq and the process repeated until the correction vector is 

small. A complete description of this elementary adjustment pro- 

40 24 

cedure is given by Richardus and Uotila , among others. 


4.22 Case (b) 

This case was first described 
sequently fully detailed in 1964. 
with the same unknown parameters, 
written as follows: 


An 42 

by Brown * in 1959 and sub- 
There 'are two observation sets 
These sets can be symbolically 


Fd/^X ) = 0 
a a 

g(l 2 ,x ) = 0 

a T a 


(4.4) 


The first set of observations results from plate observations and 
uses the collinearity conditions, equation set (2.5) j as the mathe- 
matical function. The second equation set originates from obser- 
vations on the elements of exterior orientation. These observa- 
tions are made to conform to the following concept under minimum 


variance : 



W M 
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Y . - X . = 0 
ax ax 

where X . is the adjusted or theoretical value of the parameter i, 
ai 

Y . is the adjusted observation on the parameter i, 
ax 


and i is the parameter x , y , u , ^ 


The solution of such a set of equations follows a well-defined 
routine, commencing with the linearization of the two symbolic 
matrix equations : 


E E ' 

‘A + V e i = 0 

E E 

A 2 V 2 + B 2 A + e 2 = ° 


(4.5) 


It is noted that by 
= -I. Hence the 


definition A^ = I and = I, and that 
observations can be represented by 


E E 

^ = 0 - 2n equations 

E (4.6) 

V 2 -A + € 2 = 0 - 14 equations 


The minimum variance solution of these two equations is obtained 

E 

by minimizing the variables V^, and A in the following function : 


4 


T T 

V P V +. V P V 
1 1 1 2 2 2 


m E E 

2A^(V 1 + B 1 A + € 1 ) 


- 2A 2 (v 2 - A +€ 2 ) = 0 


(4.7) 


Differentiation with respect to the variables in (4.7) yields: 
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* • i 

= 2P i v i " 2 ^l 
« 

- 2P V - 2 A 
^ ~ 2F 2 V 2 2A 2 

1 

lb ? ?Tl v 

nr = - 2 B i A i + 2A 2 

d A. 

Hence, the necessary information for a solution is as follows: 
E E 

V 1 + B 1 A + € 1 = 0 


V 2 - A + € 2 = 0 


P 1 V 1 - A 1 = 0 


(4.8) 


P 2 V 2 “ A 2 = 0 


• B 1 A 1 - A 2 = ° 


This equation set can be simplified by substituting the first 
lines into the last line from which A^ and A 2 have been removed by 
substitution of lines 3 and 4. That is, 

E E 

V 1 = - (B 1* + € 1> 

E 

V 2 = -(-A + € 2 ) 
are substituted into 

B. P„ V, . - = 0 
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which yields 


E E E E 

B 1 P l ( " B 1 A + P 2 ( “ A + Sg) 


= 0 


Simplification yields 


E m E E E„ 

(B i\ B l + V i * B i P i e i- P 2 6 2 = ° 

E E m E -i E m 

or A = -(B P x B ♦ P 2 > (B P 1 « 1 - P 2 « 2 > 


(4.9) 


The corrections are applied and a new iteration commenced with 

the updated approximate value vector, X oq . 
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Uotila has suggested the following alternative computational 


algorithm : 


Lines 2 and 3 of (4.8) may be expressed as 


A = V 2 +£ 2 


V 1 ~ P 1 ^1 


Substitution of 'these into line 1 yields 


P ^\ * V V 2 + «2> + € x = 0 


E 


(4.10) 


Similarly, line 5 of (4.8) may be expressed .as 
E„ 


a 2 - 


which can be substituted into line 4 to yield 


P 2 V 2 - k K - 


o 


(4.11) 



Regrouping equations (4.10) and (4.11) yields 


[V 1 

1 

E - 

B 1 


H 

-< 

i 


p E 

- B I«2'- 6 1 

e t 




= 


B 

-P. 


V- 


0 


2 

« 


2 




(4.12) 


The solution to this equation set in terms of V 2 is as follows: 
r E —1 ®*P I T E 1 

V 2 ' [- ( - P 2 - * 1 - P 1 V B 1 P lj [- (B 1 £ 2 ♦ *1>] 


■E„ 


E 


E 


„ ^ E 

= -CB B P x Bl ♦ Pj,)- 1 P x ( Bl e 2 + e x ) 


(4.13) 

This is normally further simplified by assuming that ^ 00 » the 
approximate value vector, equals L^, the observation vector. 
Under such assumptions, equation (4.13) reduces to 

E 


v 2 •- - (i ? p i b i + v 1 k p i e i 


(4.14) 


A new approximate value vector may be computed as + V^* 

the respective partials re-evaluated, and a new residual vector for 
V 2 computed according to (4.13). This process is continued until 
the residuals meet prescribed limits or become constant. 


4.23 Case (c) 

Case (c) is a further generalization of case (b) . This gener- 
alization is achieved by adding a third set of observations on the 
object space points. Thus the system now becomes 
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■fCl 1 , x ) = o 

a a 

g(l 2 , x ) = 0 (^.15) 

a a 1 

h(l 3 ,x ) = 0 

a a 

The last set of observations, which may be viewed as a further 
set of conditions, is extremely powerful, for now the object space 
point, the lens point and the image point are free to move so that 
the collinearity condition is fitted in a minimum variance manner 
at three points- on the ray under consideration. Since only two 
points in space- are necessary to- define a line, any additional 
points along the line do not contribute further information. How- 
ever all observations are subject to errors, and some observations 
are more easily and more precisely obtained than others. The in- 
corporation of the additional information, by minimum variance 
techniques, into a general model, allows the least precisely ob- 
served quantities to be determined with a precision approaching 
that of the quantities defining the line. Since the most diffi- 
cult and expensive observations are usually those associated with 
the object space, it is possible to relax the precision require- 
ments on these quantities by enforcing the collinearity condition 


in this general manner 
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Linearization of equation set (4.15) yields: 

EE S S 

A-jV^ + Ba + BA +£ 1 = 0 

E 

A 2 V 2 + B 2 * +€ 2 = 0 ( 4 . 16 ) 

S 

A,V_ + BA + € 2 = 0 

5 3 3 3 


which is readily simplified, by virtue of 

A x = I, A 2 = X, A ? « I, B 2 » -I, B 3 = -I, 

EE S S 

to V +3A + Ba+£ =0 

4 . 

E 

V 2 - A . + £ 2 = 0 ( 4 . 1 ?) 

S 

V z - A + €~ = 0 

3 3 



( 4 . 18 ) 

It is evident that as 0, then case (c) is reduced to 

case (b), and that as P 2 — » 0 further restrictions reduce case (b) 
to case (a). 
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4.24 Numerical Results ■ 

Since the commonly-used collinearity condition techniques may 
be applied to the numerical solution of the theoretical models 
discussed in Section 4.23, it is not considered pertinent to pre- 
sent here a detailed discussion of the numerical methods. However, 
unique features of these numerical solutions, e.g. the analytical 
differentiation of the partials, may be found in Appendix F; 

4 .241 Case (a) 

Two principal classes of tests were performed in this case. 

The first class .involved points where the exposure epoch differed 

by a constant amount from the tracking epoch. It was assumed that 

all points were exposed simultaneously by a camera system using a 

between-the-lens type shutter system. The second class of tests 

involved points imaged at differing times corresponding to the 

passage of a focal plane shutter across the image area. 

It was quickly recognised that while solutions can be obtained 

with the P = I unit weight concept, the quality of -these solutions 

leaves much to be desired. In a class one test with unit weight, 

an acceptable determination for the unknown parameters was obtained. 

However the variance-covariance matrix of the adjusted parameters ' 

24 |! -l| 

had exceedingly large norms, of the order of 10 for ||N | ^ 

Furthermore, the symmetric characteristics of this matrix had been 
destroyed. The large value for the infinity norm can be directly 
attributed to the fact that no row-column normalization was per- 
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formed prior to inversion. When the correct weight matrix, 

12 

P - 10 * I, was used, the order of the infinity norm was reduced 

12 

to 10 , while the symmetric characteristics of the variance- 

covariance matrix were retained. By applying the correct weight 
matrix, some degree of normalization is achieved, and this results 
in better computational stability. 

Class two tests all yielded very satisfactory determinations 
of the unknown parameters. The quality of the determination was, 
however, dependent on a number of factors. Of these, the most im- 
portant was the knowledge of At for a point. As indicated pre- 
viously, the image space was sectioned into eleven regions, each 
of At = 0.01. This allowed an image smear of approximately 5 
micrometers in x and 1 micrometer in y. However these values were 
not considered unreasonable, since blind velocity was a nominal 
1000 mm/sec, compared with 1200 mm/sec for Lunar Orbiter IV. The 
quality of the determination increased as At for a point became 
more refined. This refinement was achieved by reducing the slit 
width. It therefore seems that increasing the exposure interval 
by increasing the slit width is less desirable than reducing the 
velocity of the blind, even though this will increase geometric 
distortion. It was also found that the initial approximations to 
the unknown parameters had to be reasonable, otherwise no conver- 
gence occurred due to the non-linearity of the model. 

The principal results of case (a) tests are summarized in 


Table 6 



Table 6 

Summary of Case (a) Numerical Results 



Parameter 

Simulated 

Values 

Class 1 

Class 2 
smear) 
*> 

Principal 

x (mm) 

0.2000 

0.1999 

0.1996 ' 

Point 

y Q (mm) 

0.2000 

0.1998 

0.1997 

Photostation 

X (m) 

4320000. 

4320000. 

4320000 . 

.Position 

T (m) 

, 

-20p000 » 

-200000 . 

-200000. 


Z (m) 

92 

2000. 

922000 . 

921999. 

Camera 

(deg) 

0 . 

0 

0 . 

0. 

Attitude 

(deg) 

80 

.0 

80 . • 

80 . 


10 (deg) 

10 

.0 

.10. 

10. 

Photostation 

V (m/sec) 
X 

100.0 

159 . 

97. 

Velocities 

V (m/sec) 

y 

-100.0 

- 99 . 

- 106 . 


V (m/sec) 
z 

2000. 

2003 . 

2008 . 

Camera 

# (°/sec) 

0 . 

0 

O'. 

0 . 

Attitude 

£ (°/sec) 

0 . 

0 

0. 

0 . 

Velocities 

u» (°/sec) 

0 . 

0 

0 . 

.0. 

ilr 1 ! 




14 

-*#10 

— 10 13 

» BOO 








4,242 Case (b) 
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This class of tests was concerned with the effects of weighting 
the unknown parameters by observations on the parameters. It was 
possible to test with a wide range of weights, both reasonable and 
unreasonable . 

-1 


£, . * 


It is known from elementary adjustment theory that P - m 

c “to 

2 

and under the assumption that m = 1 and that 

7t ds diagonal, then P is also diagonal with elements according 
b i i ' . 

tb p u = sr = « • 

n i 

The use of reasonable weights greatly improved the stability 
of the system, when compared to case (a) . The results of an un- 
reasonable case are tabulated in Table 7, from which it is 
readily seen that excellent agreement between the simulated values 
and the adjusted values is possible. In the tabulated example, the 
position of the principal point and the attitude of the camera can 
be considered as known quantities to which almost no correction 
can be applied. It is most- enlightening to compare the norms for 
these severe cases 'with those for case (a) tests. However, due to 
the number of variables involved, it is too difficult to display 
the reduction in size as a function of increased knowledge of the 
unknown parameters. It therefore must suffice to say that drama- 
tic changes in stability occur when observations on some of the 
unknown parameters can be incorporated into the model. 
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Table 7 

An Example of Case (b) Results 



Parameter ; 

Precision 
of Obser- • 
vation t 

Weight 

Simulated 

Values 

Adjusted 

Values 

Plate 

x (mm) 
P 

io”- 5 

10 6 



Coordinates 

y (mm) 

x' 

10 -5 

io 6 : 



Principal 

x (mm) 
o 

0.001 

10 6 

0.2000 

0.2000 

Point 

y Q (mm) 

0.001 

io 6 I 

0.2000 

0.2000 

Photostation 

X (m) 

3 

io " 1 

4320000 . 

4320000 . 

Position 

Y (m) ; 

3 

io - 1 

-200000. 

-200000. 



Z (m) 

3 

io - 1 ; 

' 922000. 

922000 . 

Camera 

.ft (rad.) 

0.001 

io 6 

0.0' 

0.0 

Attitude 

<f) (rad) 

0.001 

IO 6 

80.0 

80.0 



W (rad) • 

0.001 

10? ; 

10.0 

10.0 

Photostation 

V x (m/sec) 

10 

IO -2 

100. 

100. 

Velocities 

Vy (m/sec) • 

10 

io " 2 : 

-100. 

-100 . 

r 


V (m/sec) 
z 

10 

io “ 2 ; 

2000. 

2000. 

Camera 

ft (rad/sep )j 

0.01 

4 

10 

0.0 

' 

0.0 

Attitude 

♦ 

<{> (rad/sec) 

0.01 

io 4 ■ 

0.0 

0.0 

Velocities 

•it (rad/sec)< 

0.01 

4 

10 

0.0 

0.0 


■ 

ao 




, 

10 7 
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4.243 Case (c) 

Tests oh this case were done in conjunction with the multiple 
photo tests, since n = 1 is a special case of the general n x m 
photo block. The principal’ difference is that observations are 
now available on all quantities. As would be expected from the 
previous discussion, good solutions with stable inverses are pos- 
sible under a wide range of conditions. In particular, if the 

weights of the observed quantities are correctly entered then, 

2 

since the expected value of m Q is 1, the variance-covariance matrix 

of the adjusted quantities is approximately known. The upper bound 

of the matrix norm || A b® estimated to be the quantity 

n.max(a .■}), ' since correlation is <l|l| . Consequently, for this 

i 

case, the quality of the solution can be estimated prior to execu- 
tion and, if warranted, precautions taken to ensure that an accep- 
table solution is obtained. 

The reader is referred to Section 4.3 for numerical data. 

4*3 The n x m Photo Block - The Intersection Problem 

The case of the intersection problem is built upon the single 

photo resection problem. Associated with each photograph is a set 

of observations on plate coordinates symbolically expressed as 

F(L la , X ) = 0 
a a 

F(L lb , X ) = 0 
a a 

f(l 1:l , x ) = o 

a a 


(4.19) 
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which can be expressed as 

D(L X , X ) = 0 
a a 

The same concept can be extended to observations associated with 

the elements of exterior orientation: 

G(L 2a , X ) a 0 
a a 


G(L , X ) = 0 
a a 


(4.20) 


G(L 2 ^ , X‘ ) = 0 
a a 

Again, this system can be simplified by using the matrix equation 

E(L 2 , X ) = O' 
a a 

Finally, the observations on the survey coordinates can be 
expressed as 

H(L 5 , X ) = 0 (4.21) 

a a 

• Thus, the matrix system associated with the intersection prob- 
lem may be expressed as 

DCL 1 , X ) = 0 
a a 

E(L 2 , X ) a 0 (4.22) 

a a 

H(L 3 , X ') = 0 
a a 

which corresponds to equation set (4.15) of Section 4.25. The 
solution to equation set (4.15), namely (4.l8), is. therefore also 
the solution to (4.22). However the sub-blocks are now built from 
data associated with more than one photograph. It is therefore 



74 


considered pertinent to restate the solution to equation set 
(4.19) and then to discuss the composition of the sub-blocks in 
this new context. 

ir 1 r 

E„ 


e t e i 

, B i P- L B + P ; 

1 

e t 

B P 1 

S 

B 

S m E ! 

S m ' 

S . 

m I 

B P- B ' 

T 

B P. 

B + P, 

1 i 

1 

3 

e t 

E 


The segment B 

P 1 B + 

E 2 P 


E 

A 


s 

A 


B T P 1<£l - P 2 S 2 


B P 1 6 1 - P 5 € 3 


= 0 


(4.23) 


(4.19) ’with equation set (4.20). It is noted that the elements of 
exterior orientation appear only in functions concerned with a par- 
ticular photo. That is, the differentials of elements of photo i 
are zero except in photo i. The following matrix describes the 
situation for a two photo case: 


E, 


T 


E 


B P x B + P 2 = 


B. 


’T 


1 + P- 

la 2a 


la 1 - 

,(l4x2n) (2nx2n) (2nxl4) (l4xl4) 


0 


0 


B. 


T 


E 

B, 


lb 1 lb 2b 

(l4x2n) (2nx2n) (2nxl4) (l4xl4) 


The sub-element B B + is a square 3 n x 3n matrix formed 
by augmenting equation ' set. (4.19) with (4.21). Unlike the exteri- 
or elements, ground points are not confined to imaging on a single 
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photograph. Hence, under the assumption that the points are un- 
correlated, the sub-element may be expressed as follows: 


m •"> 

B . F B + P T 

1 5 

(3nx2n) (2nx2n) (2nx3n) (3nx3n) 


B 


T 


b: 


la 2a 


_ 


' S ‘ 

P- 0 • « • 

la 


B la 



S 

0 P lb ’ * * 

• • 


B 

2a 

• 

• • 

• • 

■ 


• 

■ 


+ P. = 


S T S s s 

B n P, B + BT. P_. B , + 

la la la lb lb lb 


E T S S T E 

The sub-element B B and its transpose B F^ B combine parts 

of the previous two sub-elements. It is readily seen that the 

e t a . 

dimensions of B -P B are 

(no. of photos times 14 x 2n) . (2n x 2n) * (2n x 3n) 
and hence 

r E m S 


Em S' S„ E ~ 

B P B = (B P B) = 


5 1 

p„ 

B, 

la 

1 

la 

!* 


S 

p 

B 

lb 

1 

lb 


• 



• 



Either by expansion or by applying the above principles to the 
U section of the normal equation matrix, it may readily be deduced 
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that the following equations hold: 

®T 


B T p 1 £ 1 - p 2 e 2 = 


B i p i 6 , 
la la la 


B^. K.. € - . 
lb lb lb 


P € 

2a 2a 


P € 

2b 2b 


S S„ S ■ 

B T P_ € - P_€ = b!t P_ £ _ + B, . P. • P, € , 

11 3 3 la la la lb lb lb 3 3 


The structure of both the normal equation matrix and the con- 
stant vector are unique under the special conditions of diagonal 
weight matrices. As- the size of the matrices increases, this 
unique structure must be exploited more and more in order that the 
size of the problem does not outstrip the capabilities of the com- 
puting facilities.' In the numerical work associated with the 
testing of this model, sufficient core was available so that par- 
titioning schemes were not necessary. 


b.31 Numerical Tests , 

The numerical tests associated with this section were per- 
formed on a block of two photographs. The photographic data were 
generated according to the method described in Section k.l for 
positions and attitudes approximating Lunar Orbiter IV photographs 
102 and 109* The ground data and exterior orientation data neces- 
•sary to accomplish the data generation process were then assumed 
to be equivalent to the required observation parameters. 
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The two photographs had a nominal 60% overlap with a total of 
2k ground points appearing in the overlap area. For convenience, 
it was .assumed that the variance-covariance matrices for the obr- 
served quantities were diagonal in nature. Table 8 gives a typi- 
cal set of standard deviations associated with the observed quan- 
tities . 



Station 


X (m) 
Y (m) 
1 (m) 


3 


8 

(Continued) 
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. — 1 

' 

Observed 

Parameter 

Standard 
Deviation (£) 

Camera 

(rad=) 

0.001 

Attitude. 

(rad) 

0.001 - 


to (rad) 

0.001 

Camera 

(m/sec) 

10 

Velocity 

V (m/sec) 

•7 

10 

• 

(Vsec) 

10 

Camera Attitude 

If (ra^/sec) 

0.01 

Eates 

<j> (ract/sec) 

0.01 


w (rad/sec) 

‘0.01 





In the example given in Table 8, the rotation elements are 
considered very well-known. However, tests run with these values 
relaxed yielded similar results. It should be noted that camera 
attitude in the Apollo J missions will be well-known from the 
coupled stellar camera. 

As mentioned earlier',' the most striking feature of these tests 
was the form of .the- resulting variance-covariance matrix of the -ad- 
justed quantities. - The variances of the adjusted quantities were 
very similar to those used in the weighting matrices. This corres- 
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pondence between the variance-covariance matrix of the adjusted 
quantities and that associated with the observed quantities allows 
stability computations or estimates to be made prior to the solu- 
tion of- the system. The correct estimation of the precision of 
the observational process is therefore of great importance,. 

Consider the norm theory of Chapter 3 and in particular 
equation ( 3 * 18 ) : 


* X L - 


,-i 


03 


XI + || A 

oo 


-1 


SY 


CO 


1' - 


,-1 


oo 


E 


(3.18) 


-lb 

It is noted from Chapter 3 that IjEU^ = n.10 and that 


SY 


= 10“ 15 for the IECC 360/73' in the double precision mode. 


oo 


Then A X = N ^ and an upper limit for 

co oo 


N 


-1 


03 


can be 


estimated as n.maxCfff. ) where 6V. is ‘the variance of-an obser- 

£ 

vation. 

Similarly, NL ■ 
becomes: , 


U|l a max(u. ), hence the above equation 
co ^ i 


l Sx L 

upper limit 


n.max( ) .n.10 ^.max(u.) + n.maxC 6 ?. ) .10 ^ 
i i- i 


1 - n.max( <?? . ) .n.10 
ii 
i 


-15 


which 'yields: 


SX 


n 2 .10"’ 1 ^.max( CT. ) ,max(u. ) + n .10- ,max( Cf . ) 


ii i‘ 

i i 


ii 


upper limit 


1 - n 2 .10 -15 .max( ) 


ii 


i 



8o 


■In the tests associated with this chapter, the following condi- 
tions applied: 

n = 100 

max( O'. ) s 10^ 
i 

8 

max(u.) = 10 
x 
x 


Hence, 


||sx|L 

upper limit 


lo^.io-^ao^io 8 ^ io 2 .io~ 15 .io 6 
l - io\io- 15 .io 6 


i° 3 + io-? „ 10 3 


x.e. 


SX 


CO 


1-10 
SSr 105 


-5 


upper limit 


This value applies to corrections to the survey coordinates 
and is of the . same order as their estimated standard deviations, 
indicating that the 'adjustment procedure may not further improve 
knqwn values. 

An analysis of Jjsxlj^ after adjustment reduced its value to 

|)sx|j^ 10, which indicates the level of significance of the 

6 

solution to be approximately 1 part in 10 ... It should also be 
noted that row, column normalization was not performed on the nor- 
mal equation • matrix prior to inversion. This normalization should 
be done if the .most representative norms are to be obtained. 

In general, the procedure adopted required three iterations to 
converge to stable answers. The time for each iteration was 2 
minutes 10 seconds for a total of twenty-four points. The greater 
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part of' the time. was spent in performing the analytical differen- 
tiation required for the partial differentials. One test' was’ per- 
formed to compare ■ the time required to determine the partial by 
.numerical techniques. The following principles were applied: 


Assume that 



w a s required. . 


Then 


Hi 

bn 


(F ) 

V 1 ; 0 +S 


- (F) 


0 


where (F^)q is F^ evaluated for assumed values of the param- ■ 
eters . 

(F^)^ is F^ evaluated for the assumed values of the param- 
eters with a small delta increment added to the differential 
parameter under investigation. 

The numerical procedure was found to be an’ order of magnitude 
faster than the analytical differentiation and yielded the same 
numerical values. However the numerical agreement is dependent on 
the S increment chosen. This concept requires further investiga- 
tion as it was not feasible to continue this investigation.* 

In these block tests it was also noted that the coordinate 
values of the principal point did not seem to be as responsive to . 
adjustment as they were in the single photo tests. This was in 

part overcome by using better estimates than had been initially 

-2 

contemplated. That is, the variance was decreased from 10 to 

-3 

2 .5 x 10 for observations on these parameters. 
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In the numerical tests, it was determined that a practical 
limit of 3 iterations was necessary for the moderately perturbed 
■test data. As the perturbations became large, it was- necessary 
to complete more iterations to achieve the same level of precision. 

This level of precision, theoretically speaking, should have 
been extremely high, as it was assumed that E =1 for these 

- * ' 9. D 

tests. Hence, the residuals, V, are defined as V s L, - L =0. 

b a 

However, because of a number of factors, this will only be reached 
in the limit. The principal factors affecting this condition are 
as follows : 

(a) The augmented -collinearity -condition equations are non- 
linear - in nature. The susceptibility of these equations to ’this 
non-linearity was indicated in Section 4.241. 

(b) The weights associated with the survey stations did not 
truly represent the situation since no random errors were im- 
pressed' upon these values. Thus, the convergence rate was slowed 
down due to this incorrect weight. Ideally, weights should accu- 
rately reflect the observational precision. 

(c) The Y coordinate of the survey coordinate data appeared to 
lag behind both the X and % coordinates in reaching the 10 meter 
residual level. 

(d) The stability level of the solution was approximately 1 

6 

part in 10 , hence unit accuracy in the adjusted survey positions 
is all that may be obtained. 

Table 9 lists some typical values for tests associated with 
this section of the work. Since the precision of the adjusted - 
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Table 9 


Typical Values Associated -with the. n x m Photo Block 
with Fully Observed Parameters 



Parameter 

True 

Value 

Approximate 

Value 

Corrections I 
from first 
iteration 


x (mm) 
0 

0.200 

0.200 ! 

0.0095 


• y o (aim) 

0.225 

0.240 

-0.0149 


(m) 

0 

4320000. 

, 

4319900. 

108.1 


Y o <m) 

-200000 . 

-199990. 

-19.8. 


Z (m) 
0 

922000. 

'92190.0. 

1 

43.4 

Parameters 

W (deg) 

10. 

9-9 

■ 

o-.io 

Associated 

£ : (deg) 

80. 

79.9 

0.10 

with 

ft (deg) 

0. 

0.0 

0* 

Photostation 

V (m/sec) 

X . 

100. 

90. 

10.0 

1 

V (m/sec) 

-100. 

-95. . • 

-5.0 


V (m/sec) 
z 

2000. 

1995.. • 

5.0 


(°/sec) 

0. 

0. 

0.0 


•}> (°/sec) 

0. 

0. 

0.0 


ft (°/sec) 

0. - 

0. 

0.0 

Parameters 

x (mra) 
0 

0.200 

0.200 

0.0045 

Associated 

y (mm) 
0 

0.225 

0.240 

-0.0085 

. 

with 

X (m) 
' o 

4320000 . 

4319900. 

114.6 

Photostation 

Y (m) 
0 

200000. 

199990. 

-1.5 

2 

Z (m) 
o 

922000. 

921900. 

69.9 











Table 9 
(Continued) 
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Parameter 

True 

Value 

Approximate 

Value 


(deg) 

10. 

9.9 


(deg) 

80. 

79.9 

Parameters 

X (deg) 

0 . 

0 . 

Associated 

V (m/sec) 

X 

100. 

90. 

with 

Vy (m/sec) 

100. 

95. 

Photostation 

(m/sec) 

2000. 

1995. 

2 

w (°/sec) 

0 . 

0. 


$ (°/sec) 

0 . 

0 . 


ie (°/sec) 

0 . 

0 . 

Survey 

X (m) 

1731858.7 

1731800. 

.Station 1 

• Y (n) ; 

157282.2 

157200. 


Z (m) 

52983.5 

52900. 

Survey 

X (m) 

1696658 .2 

1696600.. 

Station 2 

' Y (m) ' 

129022.7 

129000. 


Z (m) 

344260.9 

344200. 

Survey 

X (m) 

I726572.2 

1726500. 

Station 3 

Y (m) 

-106750.7 

-106700. 


Z (m) 

192865.9 

192800. 


Corrections 
from first 
iteration 

0.10 

0.10 

0.0 

10.0 

5-0 

5*0 

0.0 

0.0 

0.0 


61.6 

46.9 

29.4 

40.2 
0.3 

63.4 

58.3 
-42.9 

44.5 
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values af£er 3 iterations has already been 'mentioned, it was de- 


cided to present instead the corrections to be applied to the 
approximate values on completion of the first iteration. Thus 
the table illustrates, in a limited manner, the convergence 
characteristics of the solution. 



5 • SUMMARY 


5.1 Conclusions 

5.11 Principal Conclusions 

The principal aim of the research presented in this report 
was to obtain, if possible, answers to those questions which were 
presented in the introductory section, Chapter 1. 

However, before elaborating on the conclusions reached by the 
adoption of the model presented in Chapter 2, it is appropriate to 
mention briefly the pitfalls encountered in two abortive attempts 
to resolve the problems. 

Preliminary attempts to correct for image motion by using 

16 

Kawachi's formulae were abandoned, as the expressions became 
unmanageable if the cause of image motion was not known. Work 
then continued with the aim of reducing this unmanageability by 
combining the individual corrections into a single uniform model, 
the existence of which had been intimated by Kawachi . "Unfortu- 
nately, this model could not be constructed. Consequently, the 
correction of image positions to a single uniform epoch correspond- 
ing to a central projection failed. 

The second set of experiments involved some recently-published 

kk 

monomorphic relationships by Das • These models included compensa- 
tory terms for focal plane shutters and image motion. However, a 
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large number of conditions exist under which the equations become 
singular. Testing of these equations indicated that these con- 
ditions would need to be completely understood before they could 
be successfully developed further. 

In general, the collinearity condition equations do not exhib- 
it singularity and, since they relate the object space to the 
image space, it was decided to augment these equations to handle 
a moving platform. The augmentation of the collinearity conditions 
with velocity and epoch of exposure does not alter the definition 
of the principal point which can be expressed as 

... the. point in a photograph or camera focal plane which 
is chosen as the centre of the image for relating the geom- 
etry of the image to the geometry of the object space. If 
the camera is distortion-free so that the geometry of the 
image is the same as that of a perspective projection of 
the object, then the principal point is the foot of the per- 
pendicular to the image plane from the centre of projection. 

k5 

National Mapping Council of Australia 
Consequently, the position of the principal point remains as des- 
cribed, since the perspective geometry remains unaltered. 

Image motion and image motion compensation do not alter the 
calibration of the interior orientation elements of the camera, of 
which the principal point is a component . However the resulting 
imagery may be deformed, such that the mathematical relationships 
between the image space and the object space are destroyed. In 
those instances when the product of exposure interval and image 
velocity is such that the expected blur is well below tolerable 
limits, the regular collinearity conditions can be used. However, 
for focal plane shutter systems where the travel time of the shut- 
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ter is considerably greater than the exposure interval, the aug- 
mented collinearity conditions should be used. In the event of the 

image motion being such that IMC needs to be applied, it is pos- 
sible to correct the observed image coordinates. This correction 
is obtained by multiplying the exposure epoch, At, by the rate at 
which IMC was applied, then summing over the epochs. 

With reference to the problem of incomplete calibration of the 
photographic system, it seems reasonable to expect that y , the 
unknown coordinate of the principal point, can be recovered dynam- 
ically from the block adjustment of the photographs. Recovery of 
this parameter is possible due to the unique nature of both the 
lunar surface and the photocoverage, resulting in very large vari- 
ations in all three coordinates across the model. It is also due 
in part to the mathematical model used, especially the incorpora- 
tion of observations on the elements of exterior orientation. 
However, there were also indications that the parameter recovered 
would be more exact if approximate values were first obtained by 
single photo space resection procedures. 

5,12 Minor Conclusions 

In addition to the solutions obtained for the main problem, a 
number of important features were observed during the work des- 
cribed above. They are as follows: 

(a) The matrix norm theory appears to offer a stable method 
with upper bounds for determining the expected and actual precision 
of matrix solution methods. The conventional P, N or M numbers, 
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on the other hand, have no supremum. The matrix norm theory ap- 
pears to offer "pre-inversion” insight for those situations where 
all parameters may be treated as observations; it also, appears 
feasible to extend this theory to include other models. 

(b) For large matrix systems, solution methods are faster and 
more stable than inverse methods. However, the variance-povariance 
matrix must be determined by secondary methods, such as the Monte- 
Carlo method. 

(c) Differentiation of complex analytical functions by numeri- 
cal methods is very much fas.ter than evaluation through their as- 
sociated analytical expressions. The accuracy of such a process 
is dependent on the 8 increment chosen for the evaluation process. 

5 .2 Recommendations 

5 .21 Recommendations Concerning Principal Conclusions 

The principal conclusions were drawn, in part, from numerical 
tests using fictitious data. It is therefore recommended that a 
small real data test be performed so that the validity of the pro- 
posed model and the conclusions are confirmed. It is not recom- 
mended that a full triangulation with Lunar Orbiter photography be 
attempted, since photography using recovered imagery will shortly 
become available from the Apollo J missions. However, many of the 
ideas expressed in this report' are directly applicable to the 
Apollo J missions and therefore warrant continued investigations. 
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5.22 Recommendations Concerning Minor Conclusions 

It is recognised that the conclusions presented in Section 
5.12 are based on insufficient data sets. There is, therefore, a 
need to increase the data base to ensure that the conclusions are 
justified. In addition, continued research should be conducted 
into the problems of stability and speed of execution, since the 
size of matrix systems in practical use continues to grow. 

It is therefore strongly recommended that theoretical and prac- 
tical research continue in the following areas: 

(a) Stability indicators. 

(b) Stable and efficient solution techniques. 

(c) Economic formation of the necessary partials. 

5.221 Recommendations on Stability Indicators 

In view of the apparent success of the norm theory in indicat- 
ing the stability characteristics of a fully-observed system, this 
approach should continue to receive attention so that stability 
indicators can be established for more general systems. In par- 
ticular, investigations should be made into the possibility of es- 
tablishing stability parameters for solutions obtained by non- 
inverse methods, such as the square root method of Section 5*6. 

5.222 Recommendations on Stable and Efficient Solution Techniques 
In Chapter 3* a number of techniques for the solution of large 

systems of equations were tested on the same large real system. 

Work of this nature must continue with all of the available algo- 
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rithms so that a complete understanding of the peculiarities of 
each algorithm can be obtained. 

Recent developments in the discipline of Numerical Analysis 
indicate that it may soon be possible to accomplish much of this 
analysis by algebraic techniques, as well as by computational, tech- 
niques . 

Thus, with continued theoretical and practical work, an order 
for computational algorithms based on size of system, desired sta- 
bility and economy may be established. This work must not be con- 
sidered outside the scope of the Geodetic Scientist and left to the 
more abstract Mathematician. 

5*223 Recommendations on the Economic Formation of Partials ■ 

In this report, a single CPU time test was conducted on the 
efficiency of numerically evaluating the required partials. This 
test indicated that numerical evaluation could be achieved in 1/10 
of the time required for analytical evaluation, without loss in 
accuracy.. It is therefore recommended that the well-known block 
triangulation procedures using the collinearity conditions be re- 
written using numerical evaluation of the partials, rather than 
analytical evaluations. The procedures are readily adaptable for 
different 9 increments and, using the appropriate statistical 
methods, the quality of the resulting solution may be compared 
against that obtained through an analytical evaluation of the par- 


tials 


APPENDIX A 


THE ROTATION MATRIX 

The well-known rotation matrix has many variants. However 
the following form is common in Geodetic Science: 


R = -M = R, R_ 
3 2 

*1 = 

E * 








■ 


• 


• 


- 


- 


* 


cos* 

sin* 

0 

• 

cos 4> 

0 

-sintji 


1 

0 

0 

i.e: M = 

-sin* 

cos* 

0 


• 0 

1 

0 


0 

COSH 

sin<*> 


0 

0 

1 


sin 4» 

0 

cos<> 


0 

-sinw 

COSO) 




■ 


. 




. 


, 


which yields on expansion:. 


M a 


cos 4> cos* cos wisin'* sinw sin* 

+ sinw sin$ cos* - cos® sin4 co sft 

-cos<j> sin* cosu cosit sinw cos-?* 

- sinw sin 4 sin* + cosw sin 4 sin* 


sin$ 


- sinw cos 4 


cosw cos^ 


It is now proposed to write the rotation matrix as a sum of two 
rotations + ft 

i.e. M = N = R/^, \ R / t i,% R ✓ . . \ 

(. rC + -fr^t) ( <P + tpCkt) ( W+ BAtJ 

i.e. N = R* R^; , R. R: R R. , 

« -Tr^t cp w u &t 


Proof that R/ , . v = R R, 

(<* +/J ) * /* 


Given that R = R = 
a 7r 


cos* sin* 0 
-sin - * cos -ft 0 


0 


0 
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and E |i = E to = 


cos ft'4t sin flat O 
-sin flat cos "Sit 0 
0 0 1 


then 


cos ft cos flat 
-sin ft sin "fiat 

-sin* cos fiat 
-cos "ft sin It at 

0 


cos ( ft + ftat) 
. -sin . ( ft + ^ at) 
0 


cos ft sin ft^t 
+ sin ft cos ft at 

-sin ft sin ft at 
+ cos ft cos ftat 

0 


sin ( ft + ftat) 
cos (ft + ftAt) 
0 


0 

0 

1 

0 

0 

1 


R (-ft+ftat) R ( * +^ 3 ) 


Similarly, it can be .proved that 


R <t> R <t>£>t “ R ( $ + < j> at) 

and 

Rw R «*at ~ R ( w +w>at) 


Hence, N may be -written explicitly in the form given on the 
following page, from which it is evident that the concept of N 
is analogous to that of the general rotation matrix M. 



c o s ( 4* + (j> At ) c o s (ir+'S At ) 


cos(w + «iAt)sin(^+ , jt&t) 

+sin (m +<o at )sin (4>+4>at ) cos (tt-KTr At) 



c o s («>+<> £*t ) co s ('£+ -<At ) 

-sin (w+iiat ) sin (4> +<}> At ) sin (*+ftAt ) 


-sin(w+iiot)cos(4»+^>At) 


sin (M>+wat ) sin ('k+H At ) 


~c o s («) + * At ) s in ( ,$ + ^'At ) c 6 s (•*+ -fc At ) 
• i 

sin(«>+*i>At ) cos ('K+'h'At ) 

+cos («» +w At ) sin (♦ +4>at ) sin (fr+ *At ) 

cos(w+wAt)cosC<J>+«j»At) 


-r 
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The rotation matrix can be differentiated readily either from 

the basic elemental rules or via skew symmetric matrices, . 

Computationally, the elemental method is superior, to the skew 

symmetric method. However, algebraically the skew symmetric idea 

46 

has many advantages (see Lucas ) . 

The following definitions shduld be noted for differentiation 
by skew symmetric matrices: 


dE 3 

doc 

d<* 

d* 


h • P 3 B 3 * h *3*3 


)» ' P 2 E 2 " E 2 P 2 


£ ' P 1 E 1 = £ Vl 


Hence, using skew symmetric matrices, the N matrix can be 
expressed as: 


* = * w* * ~k w + 


~ * P 3 R 3 R 2 B l + ^ B 3 P 2 B 2 B 1 + W E 3 E 2 R 1 P 1 


Table 10 lists the matrix components necessary to evaluate 
• • 

N by the skew symmetric method. The elemental expansion of N is 

given on page 97. 



Table 10 


Some Principal Components of the Matrices N and N 



cos* sinlf 0 
-sinlt cos'k 0 


0 

0 

1 ' 


L 


J 

L 

cos <t 

0 

-sin$ 


0 

1 

0 


sin 4 

0 

cos <f> 


r 


■ 

r 

1 

0 

0 



0 coauJ sin w 
0 -sin u cos w 


cos (icAt)] sin (if At) 0 
-sin ("*At) cos (*At) 0. 
■0 - 0 , • . 1 

i, . 

cos, ({>kt). 0 -sin ; (<*At) 

0 1 0 

» i , 

sin (<t>At) 0 cos (<{>At) 

1 ' 0 .0 
0 cos ("At) sin ("At) 

0 -sin (t<iot) cos («At). 


-ft sin (it At) ft'cos (itAt) 0 
-/(cos (K*t) -fisin (7cAt) 0 
O' 0.0 

-cj>sin (<{»At) 0 -<jicos (♦At) 

0 0 0 
4>cos ( <J> At ) 0 -$sin (^t) 

0 0 0 
0 -wsin (i^At) ucos («At) 
0 «ocos (*At) -i>sin ("At) 


0 10 

-1 0 0 

0 0 0 

0 0-1 

0 0 0 

10 0 

0 0 0 

0 0 1 

0-10 


vO 

o'* 




The N Matrix 


-4s in ( 4+ ) c o s Ck+ li&b ) 

-wsin(v+wat)sin (*+ifcat) 

c» cos (w+«iat ) sin (-ft+^cAt ) 

-* cos ( 4 +4 At ) sin (* + ) 

+*tcos (w+wat ) cosCfr+lf At ) 

+«*cos(w + wat)sin(4+4‘ i t ) cos (ft+riAt ) 

+fl?sin (w+tiat ) cos (*+■# At ) 
+cosin(w+wot)sin(4 + 4 A t)cos( , h+^tAt) 


+ 4 s i n ( <-+ wat)cos (4 + 4 il t )cos(ft+*At) 
-ftsin ( w+ w at ) sin ( 4+ 4 A- t ) sin (ft+ *at ) 

- 4 cOs(w+Wot)cOs (4 + 4 A t)cOs(‘ft+XAt) 
+£ cos (w+iiA t ) sin (4+4> A t ) sin (■»<■ +ftat ) 

4 sin( 4 + 4 A t) sin (rt+ftAt) 

-isin (w +<- At ) CO S (-Jc+fc At ) 

«COS (w>+u>At) COS (■fc+fsAt ) 

+ftc os (4+4 a t ) e° s (* +* at ) 

- 7 «c o s ( w+*i a t ) sin (“k+li a t ) 
-ucos(<"+«jAt)sin( 4+4 *t ) sin 0 * +* At) 

-^sin(w+wAt)cos( 4 > + 4 A t)sin( 7 t+A‘At) 
-Asin( w +wAt ) sin (4+4 A t ) cos +ftat ) 

-ftsin (u>+w at ) sin Oft+'tfAt ) 

-u sin ( to+<i At ) sin (4> +4 a t ) sin ('ft+'ft At ) 
+ 4 cos(»+wAt)cos( 4 + 4 A t)sin (fr+-*At) 
+H CO S ( W + lo At ) sin ( 4 + 4 At ) CO S (* + ^At ) 

4cos (4+4 At ) 

-wcos('«+tjat)cos(4+4 A t) 

+4sin(w+wAt)sin(4+4 A t) 

t 

- «sin ( <a+wAt ) c os (4 +4 *t ). 
-^cos ( w+wat) sin (4+ 4 ft t ) 



APPENDIX B 


TWO DECOMPOSITION EXAMPLES 

The decomposition of equation set (2.6) to the specific forms 
of Kawachi 1 ^’ 1 ^ is illustrated by the following two cases. 


Case 1 : Vertical photography, translation along flight direction 

The following conditions apply: 
fC = <£> = to = 0 N = I 

K = 4> = k> = 0 /. N = I 

V 


V = 


0 J 

0 
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Hence x 

P 


-f[Cx - 


- V 


(Z_ 


- v - 


(z 


-EL 


- Z q ) 


(X - X 


,>] 



f*V «(Z - z ) 

x p o 

(Z - Z ) 2 
P o 


f «V 



-f *V 
x 

h 


The minus sign is a result of the equations being formulated 

for the dispositive position with elevation up. 

• Vf 

cf . Kawachi: x^ « h being defined as modulus h. 

It is obvious that y^ = 0 since “ ^2 = [o 1 o], hence tiie 
numerator of the function equals zero. 


Case 2: Vertical photography, platform pitching (<j>) 
The following conditions apply: 
fr = £ = = 0 N = I 

i - Uj =0 

t o 

V = 0 


At B 0 
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Now for vertical photography 
of similar triangles. 





by virtue 


The y image velocity can be treated in a similar fashion. 



APPENDIX C 


MATRIX NORMS 


It is required to prove the following three identities of 
Section 3 * 1 . 

(a) Proof that (I + C) is non-singular if Jjcjj < 1. 

Consider the eigenvalue equation CX = Ax. 

Applying the rules stated in Chapter 3 to each side of the eigen- 


value . equation 

||cx|| 4 ( C | • ( x|| and |Xx|| = |A | * fx|| 

Hence, \ ± <-|| c || 


Now, (I + C) non-singular implies 


(I + C) = “JJ (1 +A.) ^ 0 


where is the eigenvalue of matrix C 


and if |a^ j< [| c| < 1, then must lie between -1 + £ and 1 - £ , 

i=n 

where € is a small quantity. Hence, tt (1 + / 0, which 

implies |(I + C) | ^ 0, and therefore the statement 

[ | < || c|| < 1 is true. 


Using the previously stated identities, it is now useful to 
derive an inequality for the matrix C = A - B. 

Consider |[ a|| = || A - B + B || || A - B| + j B || 

Then j a| - ||b|| < |a-b| 
or | A - B[| > ||A|| - | B | 
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Moreover, J| A 
and so 1 B 



|B - A | » | b|| 

Wl - I B 1 



(b) Proof 


|| (I + c) _1 | 4 



Let G = (X + C) -1 


Then. (I + C) • G = I or I = G - (-l)CG 


Taking norms in the prescribed manner, 
||l|| = ||G - (-1)CG| 

l > fla || - I| C-DcsJ 

l > |g|| - || cj • ]g| 

i » la|| • [l - ||c||] 


and hence 


J(I + C)" 1 ! « 


1 

1 - M 


(c) Proof III - (X + C)" 1 = — ™ — 

' " “ 1 - Uc|| 

Let G s I - (I + C)" 1 
Then (I + C) • G = G + CG = G - (-l)CG 
and (I + C) • (T - (I + C) _1 ), = (I + C) - I = C 
Hence C = G - (-l)CG 

Taking norms of both sides as prescribed. 


Hell > M 1 (~i)cg| 

» ||g|| - !|cg|| 

* H * [l - Sell] 

and hence ||G|| = ||l - (I + C)" 1 )! 


1 - C 


q ED - 


QED 


ilcl! 



APPENDIX D 


AN .EXAMPLE OF A TRIDIAGONAL FORM IN GEODETIC SCIENCE 
Consider the following di'fference network where station 0 is 
known : 


AIL AH AH^ AH,, AH, 

1 -2 3.4 ' 5 



T -1 ' . 24 

Then the solution is given by , -(A PA) APL t Uotila • 
Assuming P" = I, then 

1 0 0 0 0 
-110 0 0 
- A. = 0 r l 1 0 - 0 

0 0-110 

' o o- .o -i i 


(A^IA) which is 
tridiagonal 


1 - 1.0 0 0 

01-100 
0 01-1 0 

0 0 0 1 -1 

0 . 0 0 0 1 


2-1 O'O 0 
-1 ' 2 -1 0 0 
0-1 2-1 0 
0 0-1 2-1 
0' 0 0-1 1 
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APPENDIX E 


REFERENCED COMPUTER PROGRAMS 
(FORTRAN IV) 
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*44**4444 ********* t*********** 4*4*4* 4444****44*444444444*44444 4*444*4*#* 
C 

C THE REINFGRCEKENf METHOD OF INVERTING A MATRIX BY P. MORGAN JULY IS7G 
C 

C THIS IS A METHOD CF INVERTING LARGE MATRICES OUTSIDE THE CORE REGION. 

C THE DIRECT ACCESS FACILITY IS INVOKED TO HOLD THE NORMAL EQUATION 

C MATRIX AND ITS INVERSE. 

C REFERENCE FACCEEV AND FACDEEVA PP 173 

C THIS SUBROUTINE IS NOT UNIVERSAL WITHOUT CHANGES TO THE NECESSARY . 

C INTEGER CGNSTflNTS.THIS COULD BE CHANGED IF PARAMETER TRANSMISSION 

C VIA THE CALLING SEQUENCE HAS ACCEPTABLE TO THE USER. 

C UNIT 3 IS ASSUMED TO HAVE THE NORMALS WHILE UNIT A WILL HAVE INVERSE 

C THE SUBROUTINE HAS NO CALLING PARAMETERS AND DATA TRANSMI SS I ON , WHERE 

C NECESSARY IS UNDER A VARIABLE FORMAT, DOUBLE PRECISION IS USED. 

C ALL STATEMENTS PRCCEEDED BY A C/*/ CARD MUST BE MCDIFIEO TO SUIT JLB. 

C N IS ORDER CF MATRIX .THIS GOVERNS THE NUMBER OF RECORDS IN THE FILE 

C THERE BEING ONE RECORD PER ROW. 

C 

(.44 44 444 ************** *4 ***************4**4**** **444**44*444444*4**4**44*****4* 

C 

SUBROUTINE LINVRT 
IMPLICIT REAL *B IA-F,G-Z) 

C/*/ 

DIMEISSICN Al75UV(75),B(75URHO(75), 8S 175), HI 75) 

c /*/ the number of rows and/or columns of normals is oefineo 

N=75 

CALL SCL0K1 

C/4/ THE FILES ARE CEFINEC 

DEFINE FILE 3 I 75, 150 ,U,NN ), A 1 75, 1 5C, U ,NN ) 

C 

C LOAD FILE A KITH A UNIT MATRIX 

C 

DC .11 I -l , N 
DC 12 12=1, N 
A (I2> =O.CC 

12 CCNTINUC 
All) = UCC 

WRITE (A* I) I A I IA ) , I A= l , NT 
REAC (3 ' I> I A I IA ) , IA=1,N> 

IFCI.LT.N) FINC 13* t + l) 

W(n = A(I> 

DC 13 IA= l »N 

13 A ( I A ) =A I IA J/htll 

WRITE (3*1) IA(IA),1A=1,N) 

11 CCMINUE 

CC 5C ILO = I,N 

READ(3 , 11Q,ERR=1C0C > IV ( III) , 111= 1,N) 

IFIILO.LT.N) FINC (3'tlO+l ) 

vino =o.co 

CC 55 115 = UN 



10 ? 


RHGU15) = o-co 
DC 51 112 = 1 » N 

READ (4 « I 12, ERR= 100 11 ( A ( 1 13 } , 1 13= 1,N 1 
IFU12.LT. N't FIND (4*112+1) 

6(112) = A (1 15 } 

51 COM INUfc 

lFUlO.ht.I15> GC TO 45 
00 46 112 = i,N 
es( 112) = EU12) 

46 CONTINUE 
45 DC 52 112 = 1 , N 

RHCU15) = RHC (115) 46(112)* VU12) 

52 COM I hUE 
55 CONTINUE 

RHCKP1 = RF0U10) + l.DO 

DO 53 117 =l,h 

RHC U 17 )= RHC ( 117 )/ RhGKPl 

53 CONTINUE 

DO 60 120 =1,110 

REAB<4*I20,EKR=1002) ( A U 18 ) , 1 18= 1, N ) 

IFU20.LT. 110) FIND (4*120 + 1 ) 

DO 61 121 =1 , N 

A l 121 ) = A ( 1 2 l ) - RFO ( 1 2 1 ) *8S ( 120 ) 

61 CONTINUE 

WRITE(4*I20) 1AU18), 118=1, N) 

60 CCMIhUE 

TIHE=(RCLCK1 U. > J/60.0 
WR l T£ ( 6 ,998 ) 110, TIKE 
5C CONTINUE 
DG 14 1=1, N 

READ ( 3 ' I) (A(I15),I15=1,N) 

IF 1I.LT.N) FIND (3*1+1) 

DO 15 £15=1, h 

15 A(1 15)=AU15)*Vift) 

WRITE ( 3 * I) - (AU15), 115=1, N) 

REAC(4*I) (AC 115), 115=1, N) 

IFU.LT.N) F INC (4*1 + 1) 

DO 16 115=1, N 

16 A(U5)=AU15)/W( 115) 

WRITE(4*U (A( 115), 115=1, N) 

14 CChTINUE 

WR I Tfc ( 6 , 999 ) 

RETURN. 

ICCC WRlTE(6,t010) 

CALL EXIT 
1001 WHITE (6, 1011 ) 

CALL EXIT 
10C2 WRIT6(6,IC12) 

CALL EXIT 

098 FCRKAT (IX, "TUPLE NLH8ER *,I3»* NCW EXECUTED. TOTAL EXECUTION TIKE" 
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I ELAPSEC 1 N M INllTES IS SF10.3 ) 

999 FORMAT (lhtiSXi* NORMAL EXIF FRCH LINVRT ACCOMPLISHED • ) 

1010 FCRMAT < IX, 'PROGRAMMERS ERROR MESSAGE: ERROR IN READING DIRECT ACCE 
1SS DEVICE FOR V M A IR IX . EXECU T ION TERM INATED.********************* » 
2/«IX,l2CllH*I,/> 

1011 FORMAT ( IX, 'PROGRAMMERS ERRGR MESSAGE; ERROR IN READING OIRECT ACCE 
1SS CtVICE FCR A MATRIX USED TO FGRM B MATR IX. EXECUTION TERMINATED* 
2/» IX , L20 t IF * ) , / ) 

1012 FORMAT (IX, ’PROGRAMMERS ERROR MESSAGE: ERROR IN READING DIRECT ACCE 
1SS DEVICE FOR A MATRIX USED TO ACC CORRECTIONS TO. EXECUTION TERMl * 
2/ ,IX, •— NAIEC.*, 1L4(1H*>,/) 

RETURN 

ENC 
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c 

C* * 4* 4* *$«$#*$$$*** $***$***£$:(:$**:$*«#**:» S*************** 4**444* *«4*4444****«444 
MATRIX INVERSE EV REFINEMENT »P .MORGAN SEPTEMBER 1970 

THIS IS A CQUCLE PRECISION SUBROUTINE WITH 2 ENTRY POINTS FOR REFINING 
THE INVERSE LF A GIVEN MATRIX AND ITS INVERSE. 

THE PAIN JCL MUST PROVIOE FWC DISK WORK $PACES:UNITS L -AND 2 . THESF 
WORK SPACES SHCULC BE SUCH THAT TEE RECORD LENGTH EQUALS THE LENGTH 
CF A RCW. ALL CATA TRANSMISSION IS ACCORDING TO THE VARIABLE FORMAT 
ROLES. 

OISK UMTS 3 AND A FOLD THE GIVEN MATRIX AND ITS INVERSE ACCORDING TO 
RULES ESTAELISHLO FCR THE WORK SPACES. 

N IS ORDER OF MATRIX 

MAX IS MAXIMUM NUMBER OF ITTERAT ICNS 
ACC IS ACCURACY LEVEL 
MRET IS A RETURN MESSAGE CGDfc:'0,l,2 


$ $ * $ * t $ 4 4 4 44 44 444 4 4 444 44 44 44 4 3*4 444 4 444 44 4 4 44444444444 4 4 444 4 444 * 44 

SUBROUTINE REFINE ( N t KAX f ACCfMRETtNQP) 

THIS ENTRY POINT ALLOWS LOOPING WITHIN THE SUBROUTINE. DETAILED OUTPUT 
AFTER EACH STAGE GF THE CCMPUT AT ICNS . 

NOTE; IF NGP=1 THE PROGRAM SWITCHES TO EFFICIENT VERSION GF SUBROUTINE 
IMPLICIT RLAL *8 [A-H.O-Zl 
DIMENSION A(5),B(5(,C(5)iCI5) 

REhIND 1 
REWIND 2 
REWIND 3 
REWIND 4 
IREP = 0 

IF (NCF.EC.l) GO TC 716 
2 REWINO A 
READ I A > U 
REWIND A 
BB = 8(1) 

NCCUN T = 0 

1 NCCUNT = NCGUNT + l 
DO 11 1=1, N 
REAC 13) A 
DC 12 II =I,N 
REAC (A) B 

c<m=c.co 

DO 13 12 =1 , N 

Cim = C 4 1 1 ) +A ( 12 ) *G ( 12 ) 

13 CONTINUE 
12 CONTINUE 
REWIND A 
WRITEI2) C 
l l CONTINUE 
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REWIND 2 
REWIND 3 
REWIND 4 
WR IT E 1 6 » 70 1 ) 

/01 FORMAT 1 1H l» * PRODUCT AD*,//) 

DO 7C0 1=1, N 
WRITE (6 1 101 ) I 
READ (2) C 
70C hRUE(6,l02) C 
REWIND 2 
CG 15 1 = 1, N 
READ (2) C 
DC 16 11=1, N 
CUD = - CIIll 

16 CONTINUE 

cm = cm+i.co 

WRITtU) C 
15 CONTINUE 
REWIND l 
REWIND 2 
DO 17 1=1, N 
READ (l) C 
WRITE 1 2 1 C 
WRITE181 C 

17 CONTINUE 
REWIND 1 
REWIND 2 
REWIND 8 
WRITEI6, 705 ) 

705 FORRAT UH1,*UNIT NATRIX - AD »,//) 
DO 7C6 1*1, N 

WRITE16.101) I 
R FAD 12) C 

706 WR I TE ( 6 , 102 ) C- 
REWIND 2 

DC 21 I=L,N 
READ 14) 0 
DG 22 11=1, N 
CU1)=0.C0 
DC 26 12=1, N 
READ (2) C 
26 AU2)=CIIU 
REWIND 2. 

DC 23 12=1, N- 

23 DUl) = CUt)+8U2)*AU2) 

22 CONTINUE 
WRITE U) C 
21 CGMINUE 
REWIND l 
REWIND 2 
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REWIND 3 
REWIND 4 
WKITt(6,709) 

709 FORMAT (1HI, 'CHI-«D) ')- 
DO 710 1=1, N 
WRITE. (6, 101) I 
READ (1) C 
71C WRITE(6»102) C 
REWIND 1 
DO 3C I = L , N 
READ <l) C 
READ (4) B 
DC 31 11= 1 , N 
8(11) =8(11) +C (ID 

31 CONTINUE 
WR ITE (2 ) B 

3C CONTINUE 
REWIND 1 
REWIND 2 
REWIND 4 
DC 32 1=1, N 
RE AO (2) C 
WRITE (4 ) C 

32 CONTINUE 
REWIND l 
REWIND 2 
REWIND 4 

C A REFINEMENT FAS NCW BEEN ACCOMPLISHED 

C WRITE CUT TFC INVERSE 

WRITE (6, ICO) NCCUN7 
DC. 35 I = 1 , N 
WRITE (6, LCD I 
READ (4) B 

WR ITE ( 6 , 1 02 ) (8(11),IX=1,N) 

35 CCNTINUE 
REWIND 4 

IP ( NCCUNT.GT.MAX ) GO TG 36 
READ (4) E 
EBI = 8(1) 

REW INC 4 

IE ( DAeSIBB-Oen.LT. CABS (88*5.0-14) ) GO TO 36 
88 = 881 
GC TC 1 

36 DC 4C I =l,N 
0(1) = O.DC 

READ (8) A 
DG 41 Il=t,N 

D( I) = C( I) + CAES ( A ( ID) 

41 COST INGE 
4C CCNTINOE 
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REWIND 8 
AFAX = CllJ 
CO 42 I =l*N 

IF I AFAX. LT. CUD AKAX= C ( I ) 

42 CONTINUE 

OL 43 1=1, N 
0(1) = 0.D0 
READ U) A 
DO 44 11=1, N 

cm = c< i) + cabscai m> 

44 CONTINUE 

43 CONTINUE 
BMAX = cm 
DC 45 1 = 1, N 

IF ( BFAX.LT.CIl)) BMAX= DU) 

45 CONTINUE 

ACCUR = BFAX* AFAX*AFAX/ ( 1.0C-AMAX ) 

IF ( AFAX .GT. l.CC ) GO 10 60C • 

IF ( ACCUR » GT . ACC 1 GO TO 601 
HR I TE ( 6 , 107 ) NCGUNT, ACCUR 
FRET = 0 
RETURN 

601 HR 1 1 E ( 6 , 104 ) FAX , ACCUR , ACC 
IREP = IHEP + 1 
IFIIREP.GE.2) GO TC 603 
HR I T E ( 6 , 1 05 ) 

GO TC 2 
603 HR I TE 1 6, 1 06 ) 

FRET = 2 
RETURN 

60C WRITE (6, 103 ) 

FRET = 1 
RETURN 
C 

ENTRY REFIN IN, PAX ACC,' FRET , NOP ) 

C IHIS IS A SINGLE ITTERAT ICN REFINEMENT WITHOUT INTERMEDIATE" RESULTS. 

C ITS IS CONS I CER ABLY FORE EFF IC IENT, T IFEWISE, THAN THE GENERAL VERSION. 

716 CG 711 1=1 ,N 
READ {4 ) B 
DC 712 11= l, N 
READ (3) A 

cm )=c.dc 

DC 713 12=1, N 

713 c(m=cim+Ati2i*eii2) 

712 CCNUNLE 
REhINO 3 
WR 1 TE (2 ) C 
711 CONTINUE 
REWIND 2 
REWIND 3 
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RENIND A 
DG 721 1=1, N 
RE AD <4 1 6 
CO 722 1 1 = L » N 
D(Il)=C.DO 
READ (2) A 
DO 723 12=1, N 

723 D(ll) = CUll+A(I2)*8U2) 

722 CONTINUt 
REWIND 2 
WRITLCl) D 

721 CONTINUE 

REW INC 1 
REN INC 2 
REWIND A 
DO 730 1= t ,N 
READ (4) A 
READ ( l ) B 
DG 731 11=1, N 

731 D(Il)=A(Il)+A(il)-8(Il) 

WRITEI2) C 

73C CONTINUE 
REN INC l 
REWIND 2 
REN INO A 
CC 735 1=1, N 
READ (2) C 

735 NRITEU) C 
RENIND 2 
RENINU A 
KRE T=5 
RETURN 

ICC FORMAT ( 1HI,//, 1X,'*THE REFINED INVERSE MATR I X : i TTERATI ON NUMBER * , 
113,/,IX,47UE*>,///} 

101 FORMAT {// ,5X , ■ RCW NUMBER ', 1.4 ,/, 5X , * ',/) 

102 FORMAT (1X,5C25.16 ) 

103 FORMAT ( Ih 1 , / / , IX , • EXECUT ION TERMINATING: REFINEMENT NOT POSSIBUE, 
1TRY A NEW FIRST CUESS • , / , IX , 73 ( IH- ) ) 

1 04 FORMAT ( 1HI , // , IX , * F A l L EC TO REACF DESIRED ACCURACY IN SPECIFIED I 

1TTERAT ICNS.*,IX,57(1M-I,//,5X, ' [TTERA TIONS SPECIFIED = * » 1 3* / ,5X, *C 
2CM.PU TEC ACCURACY = • , C 12 . 3 , / , 5X , ' SPEC 1 F IED ACCURACY =•, G12.3,// ) 

105 FORMAT (// ,5X ,' REPEAT M0C6 INVOKED: I ITERATION COUNTER RESET TC ZER 
10.' ) 

106 FLRMAT (//,5X,* CONVERGENCE TOO SLCh : R E-EXAM IN S I TUATiCN * , / , 5X « EXEC 
1UTICN TERM IN A f INC.' ) 

107 FORMAT I // , 5X , ' S A T I S FAC TORY CONVERGENCE OBTAINED:*, /, 5X ,34 ( 1H* ) ' 

1,//,5X, 'NUMBER CF IT TERAT IONS MADE • , 12 , / , 5X , ’ AT TA INEO ACCURACY WA 
2S •, G12.3 ) 

ENC 
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A NON DESTRUCTIVE INVERSE SUBROUTINE ESPECIALLY SUITED TO POSTIVE 
DEFINITE MATRICES. NAXINUH SIZE OF This VERSION IS ICO. 

IMPLICIT REAL *8 <A-H t G-2) 

OIPENSICN At [ , h ) , B ( ifh)tP(lOQ) 

N=l-l 
' KI=P~1 
DO 1 J=1,I 
DC I K=l,I 

1 B ( J ,K ) =A(J,K) 

DC 5 K= 1,1 

DC 2 J= l,PI 

2 PtJ) = 6il,J+U/B(l,l) 

PIP) = l.CC /BUtU 

DC 4 L aI,N 
DO 3 J =I,PI 

3 BtL.J) = BIL + UJ + U- B tL+l, 1 )*Pt J ) 

A Btl*P) =~B t L+ 1 , 1 )*Pth> 

DC 5 J =1 , F 
5 B ( 1 , J ) * PtJ) 

RETURN 

END 



APPENDIX F 


ANALYTICAL FORMATION OF THE B MATRIX 
(a) Introduction 

The most difficult task in photogrammetric adjustments is the 
formation of the B matrix, more correctly termed the partial 
matrix of the functions with respect to the unknown parameters. 

■n cl function 
x.e. B = — 7 — 

<> parameter 


The mathematical functions under consideration, equation set (2.5)? 
are conveniently condensed to the following form: 



where TJ = N, (X - X - V-at) 

1 p o 

V = N„(X - X - V*&t) 

2 p o 

W = N_(X - X - V »At) 

3 P o 


Then 




£ parameter 


c> F, 


d parameter 



h TJ 

h parameter 

}V 

6 parameter 


- TJ 


- V 


jiW 

^parameter 

yw 

£ parameter 




for those cases where x q , y , and f can be considered known. 
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Thus, the problem reduces to the formation of 


h parameter 


)s parameter 


where the parameters are: 


camera attitudas4 


attitude rate change: 


7c , 

^ , (j> , ib 


photostation: 


X , Y , Z 
o o’ o 


and photostation velocities: V^, V V ^ 


(b) The Partials for the Translational and Linear Velocity Com - 
ponents of the Photostation 

Differentiation of U = W, (X - X - Vnt) yields: 

ip o 


N 0 At = -At*n 
X 0 


N -1 At = -At-n 1? 
j. 0 XeL 


au 

sf = N i ? 

o -1 


~ = -N. 0 At = -Aftt., 

z -1 . 



Differentiation of . V = N_(X 

2 p 
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-Aff , IK 

w { n ii * n 3i vr 

* F 2 -At*f 

d v x - V 

= *-JT 

o 

iF 

= 

-At'f , . U N 

W ^ n 12 + 11 32 W 

iF 

2 -At *f 

,w ~ W 
y 

. 

4 F_ 


(-n 2 i + n 31 -1 


(-n 22 + n ^ 2 w ) 


•>F 

1 _ -At*f , Uv 

dV z = W ^" n 13 + n 33 w 


bF Z -At*f , V, 

av = W ( “ n 23 + n 33 

Z 


(c) The Partials for the Attitude and Attitude Velocity Components 


of the Photostation- 

The derivation of the rotation terms is analogous to that of 
the translational and linear velocity components. 

It is again convenient to use the compressed notation of 
Section (a) of this appendix. That is, 


V = N(X - X - VAt) 

wj p 0 

- - *, - 
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Hence 


(X - X - VAt) 

p o 


It is also necessary, therefore, to form similar expressions for: 

LE 5iN 

^ ’ iu> ’ *' dw 

This can be greatly simplified by using skew symmetric matrices. 
The differentiation is as follows: 

_ ___ . R .R.*R. fc4 . *R, «R. . = P,N 
9 <frat w wat 3 


„ *At • Hi -Hi j_ »R, ,«R.. . 

55T = V Tie — * * ** 




= At«Q 1 *N where = R^-P .R^ 


= ^’^At* b$ * E ^At* R W* R W6t 


R * * t * P 2 * * R 4> At ’ Et ° R WAt 

R * * R ft At * P 2 * R ^ At * R K* R * * %a t * R 4) * R (f a t ’ R W * R W At 


= V* Where ^2 = VWV R L‘ R * 
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dN _ 


At’V^t'VVt’VV^At 

#K V%A-t •'Rp • E (j >ft . t * s u * E ^ &t * R wAt * R w* P 2 * R «*j* R w&t 

t-N.q, where ^ = B Lt’ E »V ] V E *ht 


^ N *o 

4w = \ * \at * R <fr* R $ At *Tw * R w At 


*R 


~ ^’^Afc*^ ’ R ^At* Rw * P l* R &At 


~ R * * E ^*t * H 4> * R 4 a t * R w * R *OAt * R W At * P 1 * R ttf At 


= N '^4 Where % = R «0At* P l‘ R « 


wOt 


si - vwyVt'V-dr • 4t - N - p i 


Hence the partials of U, V and W with respect to the attitudes 
and attitude rate changes are: 


*4U 





JK 

w 

4K 

= P *R-(X - X - VAt). 

1 p o 

to 

bk 

iw 


to 

L *k 


d?r 

m ^ 


At‘Q*N‘(X -X - VAt) 
x p o 
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Mj 


3U 

3<j> 


H 

iv 

= Q *N«(X - X - VAt) 

w 

d<J) 

2 p o 

d^> 



aw 



d^> 


At»N*Q_(x 
3 p 


X - VAt) 
o 


3u 


Al 

du> 


Job 

4V 

dw 

= N.Q^tXp - X Q - VAt) 

il 

d<« 

aw 



3w 


dw . 


At*N*P- (X 

1 p 


X - vat) 
o 


Substitution of these expressions into the following quotient 
rule formulae will yield the appropriate terms. 


^ F 1 

-f 

' 3U 

u 

3W 1 

*Jz 


-f 

iv 

V 

iw\ 

Zkr 

" ¥ 

1 ** 

" ¥ 


Ik 


¥ 

37: 

" ¥ 

3k | 

>F i 

-f j 

' 3u 

U 

Vw\ 

3F 2 


-f 

3v 

V 

1 


= ¥ j 

3H 

i 

“ ¥ 

wj 

w 


¥ 

37? 

~ ¥ 

37f 

j 

> F 1 

-f I 

' 5>U 

U 

i¥ 

1 * f 2 


-f 

' *v 

V 

4w\ 


_ ¥ 

1 M 

“ w 




¥ 

34> 

" ¥ 

H 

> F i 

-f 

f Vu 

u 

iw 

1 i*2 


-f I 

Al 

V 

i¥ 1 


“ W 

1 ^ 

“ ¥ 

3j 

1 '*Y 


¥ 1 

a<j> 

“ w 

j 

> F i 


1 m 

u 

3W 

^2 


-f I 

f Vv 

V 

3w 

3 u> 

_ 

" ¥ 


¥ 

)u 

J 

So) 


¥ j 

[ 

“ ¥ 

i w 

j 

>F 1 

-f 


U 

iw] 



-f 

f Jv 

V 

Aw. 1 


“ ¥ 

ids 

~ ¥ 

36 j 

ids 


¥ 

[ Sw 

“ ¥ 

3u> 
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(d) Partials Associated with the Survey Data 

It is also necessary to form the partials with respect to the 
ground coordinates for application in the total block situation. 
The same quotient rules that were used to determine the partials 

I 

with respect to the unknowns are used here with the same notation. 

Differentiation of U = N_ (X - X - VAt) yields 

1 p o 


= N 

» X p 1 


1 

0 

0 


= n 


11 


ilL - jj 

*Y 1 
p . 


o! 

l 

0 


= n 


12 


HL = n' 
b z i 

p 


o- 

0 

1 


= n. 


13 


Differentiation of V = 


V X p 


- X - VAt) yields 


2 


1 

0 

0 


= n 


21 


- K 

bY 2 
P ' 


0 

1 

0 


n 


22 


bv 

iZ. 


= N, 


0 

0 

1 


= n 


23 



123 


Differentiation of W = N^CX - X - V6t) yields 

3 p o J 



and hence 


^ F 1 -,f TJ -f 

bX~ = W (n ll ’ n 31 W 5 W = TT (n 21 " n 31 

P P 






IK 

n 32 W 







2 -f ( 

air = V <n 23 - n 33 




s:i< s:i< s;i< 
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